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CHAPTER I 
INTRODUCTION 


1.0 OVERVIEW 

This thesis Is based upon experimental and theoretical work In the area 
of gear mesh failure diagnostics. More specifically, passive diagnostic 
Instrumentation was Installed on a single mesh gear fatigue tester, located at 
NASA Lewis Research Center, to periodically record the gear mesh Induced vibra- 
tion from Initiation to failure. This information was analyzed using several 
existing gear diagnostic methodologies to determine If a correlation exists 
between the various prediction techniques and the observed modes of failure. 
This thesis presents the methods used and results obtained when applying the 
diagnostic techniques to the experimental data. 

The experimental data collected during this thesis consisted of discrete 
vibration signatures taken from the eleven gear sets that were run to failure. 
Two accelerometers, both with a frequency range of 0 to 10,000 Hertz, and an 
optical sensor, used for a shaft synchronous signal, were Installed on an 
existing spur gear fatigue tester. A timer was Installed so that the vibration 
data could be periodically recorded on an FM tape recorder. Of the eleven 
gear sets monitored, five failed by heavy wear and scoring, two failed by sin- 
gle tooth pitting, two failed by tooth breakage, and two failed by distributed 
pitting. 

The analytical work consisted of Investigating and applying several gear 
mesh failure predication techniques. Two of the methods Investigated were the 
FM0 and FM4 techniques developed by Stewart [1]. FM0 Is a general method used 
to detect a variety of failures, whereas FM4 Is more sensitive to a single 
tooth fault, such as single tooth fatigue cracks. Also used was the technique 
utilizing the Hilbert transform to demodulate the time signal. This technique, 
proposed by McFadden [2], Is predominately used to detect fatigue cracks early 
In their development. The remaining three techniques Investigated were the 
crest factor, sideband level ratio, and the non-harmonic to harmonic RMS level 
energy ratio. [3] 

This thesis is divided into five chapters to help separate the main areas 
of this investigation. The remaining two sections of the Introduction present 
the significance of condition monitoring in rotorcraft, and a brief overview 
of various transmission condition monitoring methods. Chapter II covers basic 
signal processing theory and techniques used in this analysis, and specific 
theory on each of the diagnostic methods investigated. Chapter III addresses 
the experimental work conducted during this study Including a detailed descrip- 
tion of the various gear failures found, and the results of frequency response 
measurements performed on the fatigue test rig. Results of applying the vari- 
ous gear diagnostic methods to the data collected are correlated to the actual 
gear failure modes in Chapter IV. Finally, Chapter V presents the conclusions 
of this research and recommendations applicable to further Investigations into 
gear diagnostics. 

1.1 SIGNIFICANCE OF GEAR CONDITION MONITORING IN ROTORCRAFT 

In aerospace applications, where weight and size are premiums, gear sys- 
tems are usually designed close to their projected limits. For rotorcraft 
transmissions, this design constraint translates into frequent transmission 
overhauls and frequent transmission related accidents resulting in death, 

Injury, costly damage, and fleet grounding. Current on-board condition moni- 
toring systems not only provide Insufficient time between warning and failure, 


but also result In many false alarms causing unnecessary and costly repairs. 
Thus the need for a reliable gear train condition monitoring system Is para- 
mount to the increased safety and cost efficient operation of current and 
future rotorcraft. 

Some specific examples can be given illustrating the effects of transmis- 
sion related problems on fleet readiness and safety. The entire Marine Corps 
fleet of 93 CH-53E Super Stallion Helicopters were grounded In February, 1987, 
due to a defective gear. [4] Also, The United Kingdom Department of Transport 
Air Accidents Investigation Branch (AAIB) released an official accident report 
on the November 1986 accident involving a Boeing 234 Commercial Chinook in the 
North Sea, which killed 45 people. [5] The main cause of the accident was 
determined to be a catastrophic failure of a modified gear In the forward 
transmission. This failure led to the desynchronization of the twin rotors 
such that the forward and aft rotors clashed. 

The early detection of incipient failure in a rotorcraft transmission is 
not only Important for safety, but also overall maintenance is Improved and 
becomes more cost efficient. Currently, rotorcraft transmissions are over- 
hauled at prescribed intervals. The transmissions are overhauled long before 
the calculated life of any of the critical components. Implementation of a 
reliable on-board condition monitoring system can Increase the time between 
overhauls, or on-condition maintenance, without sacrificing safety by continu- 
ously monitoring the state and wear condition of critical components. The 
reduction in the frequency of overhauls not only decreases maintenance costs, 
but also decreases the probability of replacing good components with defective 
ones. 

1.2 OVERVIEW OF TRANSMISSION CONDITION MONITORING METHODS 

There are currently two basic methods of monitoring the condition of drive 
train components. The first uses a debris monitoring device that detects the 
size and rate of wear particles in the transmission lubricant as an indicator 
of severe wear and incipient failure. The second method Is based upon vibra- 
■tion data obtained using one or more sensors mounted on the transmission case. 

Most debris monitoring devices use magnetically captured debris to detect 
surface fatigue failures in critical gearbox oil wetted components such as 
gears and bearings. Some of the debris monitoring devices classify the cap- 
tured particles by rough sizes, and keep account of the rate of capture and 
total debris count. This provides an indication as to the damage severity of 
the falling component(s) . One problem with these devices Is their limited 
ability to detect gear failures, since in most cases these-fai lures do.not 
produce large amounts of metallic debris far enough in advance to provide suf- 
ficient warning time. 

The use of vibration analysis methods for condition monitoring can be fur- 
ther classified into time domain and frequency domain methods. Time domain 
methods use statistical analysis techniques on direct or filtered time signals 
to detect parametric or pattern changes as transmission components wear. 
Statistical methods such as standard deviation and kurtosis are used to qualify 
general wear from tooth specific damage, respectively. Frequency domain 
methods use the Fast Fourier Transform (FFT) to convert the time signal into 
its 1 corresponding frequency components. Vibration energy at specific frequen- 
cies (i.e. primary, harmonics, sidebands, etc.) can be used to monitor gear- 
train component failures. Both methods use time synchronous averaging to 
cancel out all vibration that is non-periodic with the shaft frequency being 
used as the synchronous signal. 
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CHAPTER II 

THEORY AND TECHNIQUES 

An understanding of the theory behind the diagnostic techniques Investi- 
gated along with some basic digital signal processing theory was needed to suc- 
cessfully carry out this analysis. A review of this study Is presented In 
this chapter to give the reader sufficient background Into the methods used. 

2.0 DIGITAL SIGNAL PROCESSING TECHNIQUES 

An Important part of the experimental portion of this study was the con- 
version of the analog signals recorded Into digitized data. Digitizing the 
data allowed the techniques being investigated to be applied to the data effi- 
ciently and easily on a standard personal computer. Along with the analog to 
digital conversion, the digital Fourier transform was used to convert digitized 
time data to frequency data. Because of the significance of these and other 
digital signal processing techniques used to the overall analysis, some digital 
signal analysis theory will be presented below. ; 

The goal of the analog to digital conversion process Is to obtain this 
conversion while maintaining sufficient accuracy In frequency, magnitude, and 
phase Information. Two basic concepts responsible for the accuracy of the con- 
version are sampling and quantization. Sampling Is the rate of time signal 
digitization, or the timing between the digital pieces of the time record. 
Quantization Is related to converting the analog amplitude value to a digital 
value. Of the two concepts, sampling Is usually the more critical because It 
effects frequency as well as amplitude accuracies. 

All analog to digital (ADC) conversion devices sample at constant Incre- 
ments of time. Two theories, Shannon's Sampling Theory and Rayleigh's criteri- 
on, relate the sampling process to the limits of the Information obtainable 
with minimum errors. Shannon's sampling theory, also known as the Nyqulst cri- 
terion, states that the frequency the ADC device samples must be greater than 
two times the maximum frequency of interest. This theory outlines the highest 
frequency detectable In a signal, with minimum errors, given a specific sam- 
pling rate. Raleigh's criterion Is related to the ability to resolve closely 
spaced frequency components. For a given time record of length T (seconds), 
the minimum resolution frequency, or the lowest frequency component measurable, 
is given by the following equation. 

f(mln) - y <1> 

The errors associated with sampling processes can be grouped Into two cat- 
egories, variance and bias. Variance Is due primarily to random deviations of 
the sampled signal from the mean, or expected value. These random deviations 
are often caused by noise In the signal and measurement equipment. The easi- 
est way of reducing variance errors Is by time synchronous averaging. With 
this method the signal Is averaged in the time domain with each record start- 
ing at the same point in a cycle, as determined by a synchronous signal. The 
portions of the signal not coherent with the synchronous time period, such as 
random noise and signals not periodic within this period, are averaged out. 
Ideally the number of averages improves the signal to noise ratio by the fol- 
lowing equation. 

S/N (after Nt averages) = ( VNt )*(S/N) (original ) (2) 

where : S/N = signal to noise ratio 
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Bias errors cannot be minimized through averaging techniques. Limitations 
of the equipment or techniques used are the prime contributors of bias errors. 
Aliasing, a common bias error, results from the limitation of using finite 
Information to describe an Infinite process. An example of aliasing can be 
seen In In Figure 2.1. The solid line plot represents the actual signal. 
Because the sampling rate was conducted at a frequency lower than the frequency 
of Interest, the digital representation of the signal, the dashed line plot, 
Incorrectly Implies a low frequency signal. To guard against aliasing prob- 
lems, Shannon's sampling theory along with anti-aliasing filters, or low pass 
filters, should be utilized with any ADC device. The anti-aliasing filters 
serve to assure that any frequencies higher than the frequency of Interest are 
removed from the signal prior to the ADC process. 

One of the most powerful tools used In signal analysis Is the Fourier 
transform. The basic theory behind the Fourier transform Is that any continu- 
ous time signal can be represented by a series of sine and cosine functions of 
various amplitudes, frequencies, and phase relationships. Use of the Fourier 
transform produces no new Information, It just transforms the data from one 
Independent variable to another. In most instances time domain data Is trans- 
formed to the frequency domain. Determining the frequency content of a time 
based signal Is especially useful when analyzing mechanical systems with 
rotating components. The Fourier transform in the classical case Is given In 
the following equations. 


F(w) = f (t)exp(-jwt) dt (forward transform) 


f(t> >= F(w)exp( jwt) dt (Inverse transform) 


where : f(t)= function in the time domain 

F(w)= function in the frequency domain 

Because the limits of integration of the classical form of the Fourier 
transform are from negative to positive Infinity, the Discrete Fourier Trans- 
form (DFT) must be used with experimental data. Equations for the DFT are 
given below. 


N- 1 

F(m)= 1/N ^ f (n)exp(-j2irmn/N) (forward transform) 

n=0 

N-l 

f (n)= 1/N 2 F(m)exp(j2irmn/N) (Inverse transform) 

m=0 


The DFT Is based upon a set of assumptions concerning the sequence of the 
discrete data values collected. These assumptions are: 1) the signal must be 
a totally observed transient within the time period of observation, or 2) the 


4 


signal must be composed only of harmonics during the time period of observa- 
tion. Leakage, a bias error, Is a good example of the type of errors resulting 
if neither of the assumptions are met when applying the DFT. Figure 2.2 Illus- 
trates the resulting transform, labeled "DISCRETE FT", as compared to the true 
transform, labeled "EXACT FT". The time signal In this case was neither tran- 
sient nor periodic within the time record, thus the energy in essence leaked 
out into other frequency lines not representative of the actual signal. Win- 
dowing Is one method used to minimize leakage errors. For example, the Hanning 
window multiplies the time signal by a function that gradually reduces to zero 
at the end of the time period. This forces the signal to appear to be tran- 
sient within the time period, thus satisfying the DFT assumptions. There are 
functions that naturally meet the DFT assumptions, these are called self win- 
dowing functions. These functions produce no leakage errors and thus require 
no wi ndowl ng . 

The Fast Fourier Transform (FFT) was developed to reduce the computational 
time of the standard DFT process. The FFT time savings results from partition- 
ing the DFT matrix Into a number of smaller Individual matrices and reordering 
the data for more efficient computation. This results In the computational 
time required for the FFT to be a function of Nlog 2 N instead of Ns for the 
DFT. Along with the time saving advantages of the FFT Is the retention of the 
phase Information, which permits the evaluation of multichannel measurements 
such as transfer functions, and coherence. [6] 

The frequency response function is another Important digital signal pro- 
cessing method which Is used to describe the dynamic properties of a physical 
system. The classical definition of the frequency response function is the 
Fourier transform of the unit impulse response function, as shown in the fol- 
lowing equation. [7] 


H(f )= h(t)exp(-j2irft) dt (5) 

Jo 

where : H(f)= Frequency response function 

h(t)= Unit impulse response function 


The frequency response function gives a direct relationship of the output to 
the Input of a system In the frequency domain. This relationship is especially 
useful when trying to determine how a signal is altered as It passes through a 
system. 

There are several ways of estimating the frequency response function using 
measured quantities. The first method is simply computed by dividing the 
Fourier transform of the output by the Fourier transform of the input, as given 
in Equation 6. 

Hyx(f)= Y(f)/X(f) (6) 

where : Hyx estimates H 

Y(f ) and X(f) are the Fourier transforms 
of the output and input, respectively. 
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To reduce the variance error In the estimate of the frequency response a number 
of samples must be involved, as seen in Equation 7. 



(7) 


The above estimates are reasonable, however, to obtain the optimum estf 
mate of the frequency response function in the presence of noisy s1 9 jals, t 
least squares technique is often used. This technique produces an unblase 
estimate of the frequency response function even in the presence of Input • 
noise assuming that the input signal and the input noise are not correlated. 
[7] The frequency response estimate in this case 1$ given by Equation 8. 


Hyx(f) = Gyx(f) / Gxx(f) 


( 8 ) 


NS 

where: Gyx(f) * E Y(f) 1 X*(f) i 

1 = 1 


NS 

Gxx(f) - £>>, X*<f), 
1 = 1 


In this equation Gyx represents the cross power spectrum, and Gxx 

the auto power spectrum of the input, both of which can be calcu at 

the Dresence of noise. The phase information is preserved In the term Gyx. 

The cross and auto power spectrums are also used to calculate another Important 
signal processing parameter, the coherence function. _ 1)f rpl , 

The coherence function provides a means of determining the cau y 
tionship between the input and output of a system. The coherenc 
given In Equation 9 below. 


Y 2 (f) 


Gxv(f) Gxv*(f) 
Gxx(f) Gvv(f) 


(9) 


where: 


Gxy(f) = The cross spectrum 

Gxy*(f) = The complex conjugate of the cross spectrum 
Gxx(f) = The auto power spectrum of the Input 
Gyy(f) = The auto power spectrum of the output 


The coherence function is a measure of the amount of the output signal that is 
directly related to the input signal at any specific frequency. If , 1 for 
example , at a certain frequency the coherence function has a value of 1.0. then 
the output is a direct result of the input signal at that frequency. 
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coherence function has a value less than unity but greater than zero, then one 
or more of the following conditions may exist, as outlined by Bendat [7]: 

1) Extraneous noise Is present In the measurements. 

2) Resolution bias errors are present In the spectral estimates. 

3) The system relating y ( t > to x(t) Is not linear. 

4) The output y(t) Is due to other Inputs besides x(t). 

The coherence function Is normally used In conjunction with frequency response 
measurements. Here the coherence function provides a measure of the quality of 
the measurement by giving the relationship of the output signal to the Input 
Impulse. Another way of representing the coherence function Information Is by 
the signal to noise ratio. The signal to noise ratio Is easily found using the 
coherence function as shown In Equation 10. 

2 

S/N (f) = Y (10) 

1 - y (f ) 


2.1 THEORY OF GEAR FAILURE PREDICTION METHODS 

Several gear failure prediction methods were Investigated and applied to 
the experimental data. The basic theory behind each method Is given In the 
following sections. 

FMO METHOD 

The FMO parameter Is a time domain discriminant proposed by Stewart [1] 
that provides a simple method to detect major changes In the meshing pattern. 
FMO is defined as the ratio of the peak-to-peak level of the signal average to 
the sum of the rms levels of the meshing frequency and Its harmonics, as given 
In Equation 1 1 . 


FMO 


PP 

2 A(f.) 

1 = 1 1 


( 11 ) 


where: PP = peak-to-peak level of signal average 

A = amplitude at mesh frequency (1=1) and 
harmonics ( 1 > 1 > 

FMO Is formulated to be a robust Indicator of major faults In a gear 
mesh. Major tooth faults, such as breakage, usually results In an Increase In 
peak-to-peak level and no appreciable change In meshing frequency levels, 
which causes FMO to Increase. For heavy uniform gear wear, the peak to peak 
level remains somewhat constant and the meshing frequency levels tend to 
decrease, causing the FMO parameter to Increase. In heavy wear situations, 
the gear tooth profile and surface quality degrades, causing the meshing energy 
to be redistributed from the meshing frequencies to the modulating sidebands. 
From the physical viewpoint, two surfaces In sliding contact produce more ran- 
dom fluctuations In the signal as the amount and magnitude of their surface 
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variations Increase. This explains the redistribution of energy In the fre- 
quency domain with the increasing wear. For minor tooth damage, neither peak- 
to-peak levels or meshing frequencies levels change an appreciable amount, 
thus FMO response In this case Is limited. C81 

FM4 METHOD 

FM4 was developed to detect changes In the vibration pattern resulting 
from damage to a single tooth. FM4 analysis, proposed by Stewart [1], filters 
out the regular meshing components from the signal average and performs two 
statistical operations, standard deviation and kurtosls, on the difference sig- 
nal. Equation 12 and Figures 2.3 through 2.5 illustrate the formulation of the 
difference signal . 

D(t) = A ( t ) - R( t) (12) 

where: A(t> * original signal average 

R(t) = regular meshing components of signal average 
D(t> = difference signal 

Figure 2.3 represents a plot of an actual signal A(t) after time synchro- 
nous averaging. Figure 2.4 represents a plot of the regular meshing components 
R(t) of that signal. R(t> is found by taking the FFT of the original signal, 
extracting the regular components and taking the Inverse FFT of these compo- 
nents. The regular components consist of the shaft frequency and Its harmonics, 
the primary meshing frequency and Its harmonics along with their first order 
sidebands. The first order sidebands are generally due to runout of the gear 
because of machining or assembly Inaccuracies [9], and thus are considered 
regular components. The difference signal Is then found by subtracting the 
regular components from the original signal. Figure 2.5 shows a plot of the 
resulting difference signal. The FM4 analysis method Is comprised of the 
standard deviation and kurtosis of the difference signal along with a squared 
representation of the difference signal. The square of the difference signal 
represented In Figure 2.5 Is shown in Figure 2.6. As seen In this figure the 
square of the difference signal magnifies any abnormalities present in the dif- 
ference signal . 

The FM4 analysis method uses standard deviation and kurtosls to extract 
information from the resulting difference signal. The standard deviation of 
the difference signal Indicates the amount of energy In the non-meshing compo- 
nents, where the kurtosis indicates the presence of peaks in the difference 
signal. The theory behind FM4 is that for a gear In good condition the dif- 
ference signal would be primarily noise with a Gaussian amplitude distribution 
[81. The standard deviation should be relatively constant, and a normalized 
kurtosls value of three would be reflected [3]. When a tooth develops a major 
defect, such as root cracks or severe spalling, a peak or series of peaks 
appear In the difference signal. The kurtosis value would thus Increase to 
reflect these peaks. The standard deviation will Increase only when the 
peak(s) become severe enough to bring up the rms level of the entire difference 
signal. It should be noted that the standard deviation will also Increase with 
increases In uniform gear wear. This is to be expected, since the standard 
deviation of the difference signal is an indicator of the energy level of the 
non-meshing components, which increases as the gear wears. Thus both the 
standard deviation and the normalized kurtosis should be used in conjunction 
with each other to detect the onset and severity of single tooth defects. 
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Standard deviation Is defined as the second statistical moment of an 
array of values about the mean of those values. It Is an Indicator of the var- 
iance of the values In the array. The standard deviation of the difference 
signal can be found using the following equation. 


RMSDS = 




(13) 


where : RMSDS = Standard Deviation 

d - Mean value of the signal 


Kurtosls Is defined as the fourth statistical moment of an array of values 
about the mean of those values. It Is an Indicator of the existence of major 
peaks In the array. The digital form of the kurtosls equation Is given In 
Equation 14. 



1 = 1 

where: K = Kurtosls 

d = Mean value of signal 


(14) 


This absolute kurtosls value will Increase proportionally with general Increases 
In the standard deviation. To keep the kurtosls parameter sensitive to single 
tooth damage only, the normalized kurtosls equation Is given below. 


NK = 



(15) 


where: NK = Normalized kurtosls 

d = Mean value of signal 


As seen In Equation 15, the normalized kurtosls Is simply the absolute 
kurtosls divided by the fourth power of the standard deviation, or the square 
of the variance. This allows the normalized kurtosls to be sensitive only to 
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peaks In the difference signal, regardless of changes In the standard 
deviation. 

Because the normalized kurtosis values are non-dimensional, signals witn 
different magnitudes but similar shapes will have similar values. A square 
wave Is found to have a normalized kurtosis value of 1.0, for a sine wave the 
value is 1.5, and for a signal of essentially noise with a Gaussian amplitude 
distribution the normalized kurtosis value Is found to be 3.0 [3]. Thus, a 
normalized kurtosis value greater than 3.0 Is Indicative of a peak or series of 
peaks existing In the signal. One pitfall with the normalized kurtosis parame- 
ter Is Its drastic decrease in peak sensitivity as the number of peaks of simi- 
lar magnitudes increase beyond two. Thus for failures Involving two or more 
teeth, the normalized kurtosis value may not increase far beyond 3.0, and the 
failure must be detected with the standard deviation level. 

HILBERT TRANSFORM METHOD , , . _ 

A technique was developed by McFadden [2] to detect local gear defects, 
such as fatigue cracks, using the Hilbert transform. The basic theory behind 
this technique 1s : that the sidebands around the dominant meshing frequency 
modulate the meshing frequency to produce the time average signal. Using the 
Hilbert transform, the signal can be demodulated, resulting In the correspond- 
ing amplitude and phase modulation functions. The phase modulation function 
is especially sensitive to fatigue cracks by indicating a phase lag at the 
point the cracked tooth goes into mesh [2]. 

The Hilbert transform is primarily used to transform a real time signal 
into a complex time signal with real and imaginary parts. The real part of 
the complex time signal is the actual time signal, and the imaginary part Is 
the Hilbert transform of the actual time signal. This complex time signal is 
referred to as the analytic signal, and is given In the following equation. 


AN(t) = A<t) + i H[A(t) ] 

where: AN(t) = analytic signal 

A(t) = original signal 

HCA(t)] = Hilbert transform of original signal 


The Hilbert transform of a real valued time signal , A(t) , Is defined as the 
convolution of the time signal with 1/irt, as shown in Equation 17. 


H[A( t) 1 = 1 



A(x) ( 1 /( t-x) ) dx 


(17) 


An easier way of computing the Hilbert transform is to utilize the convolution 
theorem, i.e. the convolution of two signals in the time domain is equivalent 
to the inverse Fourier transform of the product of the two signals in the fre- 
quency domain. Thus, to determine the Hilbert transform of a real-valued time 
signal the signal must be transformed to the frequency domain, phase shifted 
by -90 degrees, and transformed back to the time domain. The determination of 
the Hilbert transform using this method is illustrated In Equation 18. 
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(18) 


H[A(t) ] - F"U(-1 sgnf) A(f)] 

where: A(f) = Fourier transform of original signal 
sgnf = 1 for f>0 , -1 for f<0 
F - U] = Inverse Fourier transform 


Once the analytic signal Is found, the amplitude and phase modulation 
functions can be determined. The amplitude modulation function, or envelope 
of the signal. Is the magnitude of the analytic signal, and Is given In 
Equation 19. 


| AN( t) | = [(A(t)) 2 + (H[A(t>]) 2 ] 1/2 (19) 


The phase modulation function, or instantaneous phase variation, Is found by 
using Equation 20 below. 

<j>(t) = tare 1 CHCA(t) J/A(t)] -2irf 0 t (20) 

where: f 0 = carrier frequency 


The phase modulation function, 4>(t), represents the Instantaneous phase angle 
variation with respect to the nominal gear speed CIO]. The second term In 
Equation 20 represents a ramp function with a frequency equal to the carrier 
frequency, f 0 , that is being modulated. This term Is required to separate 
the Instantaneous phase angle variations from the constant carrier frequency 
phase function. 

Before creating the analytic signal, the original signal must be bandpass 
filtered about a dominant meshing frequency. This dominant frequency Is either 
the primary mesh frequency or one of its harmonics, whichever appears to give 
the most robust group of sidebands. The width of the bandpass filter depends 
on the location of the meshing frequency to other meshing frequency harmonics. 
McFadden [2] suggests using a bandwidth giving the maximum amount of sidebands, 
even if they interfere with the sidebands from other harmonics. The reasoning 
here Is to Include as much of the primary modulating sidebands as possible, 
assuming that the interference from the other sidebands Is negligible. The 
dominant frequency, or carrier frequency, is then removed and the resulting 
bandpassed sidebands are used in constructing the analytic signal and modula- 
tion functions. 

As stated earlier, the phase modulation function appears to be the most 
sensitive of the two modulating functions to gear fatigue cracks. One particu- 
lar feature of the phase modulation function Is Its ability to distinguish 
between an actual fatigue crack and a particle on the tooth. For a fatigue 
crack, the phase modulation function would experience a phase lag at the point 
the cracked tooth went into mesh. For a particle on the tooth, the phase modu- 
lation function would experience a phase lead when that tooth went into mesh 
[ 2 ]. 
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CREST FACTOR . , 

The crest factor, CF, is a simple measure of detecting changes in the sig- 
nal pattern due to impulsive vibration sources, such as tooth breakage [3]. 

The crest factor Is easily calculated by dividing the peak level of the signal 
average to the standard deviation of the signal average, as given in 
Equation 21 . 



where: PL - peak level of signal 

RMS = standard deviation of signal average 


( 21 ) 


Peaks In the signal will result in an Increase in CF. 


SIDEBAND LEVEL FACTOR 

The sideband level factor, SLF, is a course Indicator of single tooth 
damage or gear shaft damage. [8] To calculate SLF the first order sideband 
levels about the primary meshing frequency are divided by the standard devia- 
tion of the signal average, as seen in Equation 22. 


SLF 


FOSL 

RMS 


where: FOSL = first order sidebands level 

RMS = standard deviation of signal 


( 22 ) 


A bent or damaged shaft will result in an eccentric meshing pattern, directly 
reflecting an Increase in the first order sidebands level, thus increasing the 
SLF value. It is unclear to the author how this factor is sensitive to single 
tooth damage. 

ENERGY RATIO 

The energy ratio (ER) is formulated to be a robust Indicator of heavy uni- 
form wear. This factor divides the standard deviation of the difference signal 
by the standard deviation of the signal composed of the regular frequency com- 
ponents only [3]. The difference signal and the regular components signal are 
the same as those defined previously in the FM4 analysis section. Equation 23 
illustrates the ER factor. 


RMSDS 
= RMSRC 


(23) 


where: RMSDS = standard deviation of difference signal 

RMSRC = standard deviation of regular components 
portion of the signal 
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For heavy wear, the mesh energy redistributes from the regular meshing fre- 
quency components to the modulating sidebands. This would result In an 
Increase In the rms level of the difference signal, the numerator, and a 
decrease In the rms level of the regular components signal, the denominator. 
The net result would be an Increase In the ER value In a more robust fashion 
than just the Increase In the difference signal rms level, as calculated In 
the FM4 analysis. 
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Figure 2.1 : Example of the Aliasing Problem. 



Frequency 


Figure 2.2 : Example of Leakage Problem. 
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Figure 2.3 : Plot of Actual Signal Average, A(t), for One Shaft 

Revolution. 



Figure 2.4 : Plot of Regular Part, R(t), of Signal, for One Shaft 

Revolution. 
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PHAPTFR TTT 

EXPERIMENTAL PROCEDURES 


The experimental portion of this study was conducted on a gear fatigue 
testing apparatus at NASA Lewis Research Center. Some background Information 
on the fatigue tester will be presented In this chapter along with some details 
on the type of gears tested and their specific modes of failure. Also dis- 
cussed In this chapter is the test instrumentation used, basic test set-up, and 
some frequency response measurements performed on the fatigue tester. 

3.0 SPUR GEAR FATIGUE APPARATUS 

The spur gear fatigue test apparatus at NASA Lewis Research Center first 
became operational in 1972. It has supplied the means to obtain crucial data 
on the effects of gear materials, gear surface treatments, and lubrication 
types and methods on the fatigue strength of aircraft quality gears. A cut 
away view and schematic diagram of the fatigue rig Is shown In Figure 3.1. The 
test rig operates on the four-square, or recirculating torque, principal. As 
seen In Figure 3.1, the slave gear provides the recirculating torque by virtue 
of the hydraulic pressure acting on the load vanes. The tooth load Is con- 
trolled by varying the amount of hydraulic pressure applied. Because of this 
recirculating torque principal, the drive motor provides only the power 
required to overcome the frictional losses at the desired running speed. The 
fatigue rig is capable of providing 75 kW (100 hp) at the designed operating 
speed of 10,000 rpm. 

3.1 TEST GEARS 

The gears tested using the fatigue test apparatus were all the same type 
and were subjected to identical loading conditions. All of the gears were made 
of AISI 9310 steel and were manufactured to AGMA class 13, aircraft quality, 
tolerances. The test gears have 28 teeth, a pitch diameter of 88.9 mm 
(3.50 In.), a pressure angle of 20 degrees, and a tooth face width of 6.35 mm 
(0.25 in.). The gears were loaded to 74.6 Nm (660 in lbs), which resulted In 
a pitchline maximum hertz stress of 1.71 GPa (248 ksl), at an operating speed 
of 10,000 rpm. This represents a load of almost two times the normal design 
load for these gears. This testing procedure was used In order to obtain fail- 
ures within a reasonable amount of test time. Table I lists the basic run 
parameters along with the number of hours to failure and the overall mode of. 
failure for each run. As seen In this table, the only run parameters that dif- 
fer between the runs are the surface treatment, and the type of lubricant. 

Three different surface treatments were encountered in the eleven runs 
recorded. The gears in runs 1 through 7 all used the standard surface treat- 
ment, which is a final grinding operation to a surface finish of 356 mm rms 
(14 in. rms). The gears in runs 8 through 11 were shot peened and then honed. 
The main reason for shot peening is to create subsurface residual stresses that 
improve the pitting fatigue life of a gear. The gears In runs 8 through 10 
used the SPH method, a high intensity shot peening method designed to produce 
high subsurface residual stresses. The gears in run 11 used the SPL method, a 
low Intensity shot peening method designed to produce low subsurface. residual 
stresses. The basic differences in life and failure modes, as seen in Table I, 
are primarily due to the different lubricants used. 
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TABLE I: Test Run Descriptions 



GEAR SURFACE 


RUN ft 

TREATMENT 

LUBRICATION 

1 

Standard 

Type A 01 1 

2 

Standard 

Type A 01 1 

3 

Standard 

Type A 01 1 

4 

Standard 

Type A 01 1 

5 

Standard 

Type A 01 1 

6 

Standard 

Type A 01 1 

7 

Standard 

Type A 01 1 

8 

SPH 

Type B 01 1 

9 

SPH 

Type B Oil 

10 

SPH 

Type B 01 1 

11 

SPL 

Type B 01 1 


OTAL LIFE 
(Hours) 

BASIC 

FAILURE MODE 

49 

Heavy Wear 
and Scoring 

36 

Heavy Wear 
and Scoring 

9 

Broke Tooth 

79 

Heavy Wear 
and Scoring 

39 

Heavy Wear 
and Scoring 

54 

Heavy Wear 
and Scoring 

50 

Broke Tooth 

520 

No Failure 
(Single Pits) 

245 

Single Pits 

339 

Distributed 

Pits 

100 

Distributed 

Pits 
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The two different oils used were classified type A oil and type B oil. 

The type A oil did not appear to have sufficient additives to provide a good 
Elasto-Hydro Dynamic (EHD) film thickness between the gear surfaces under 
extreme pressure. Without an adequate EHD film thickness the tooth surfaces 
obtain metal to metal contact, causing severe surface wear In a relatively 
short period of time. This Is especially evident by the fact that the runs 
with the type A oil, runs 1 through 7, experienced heavy wear and surface scor- 
ing, and an average life of only 63 hours. Comparatively, runs 8 through 11, 
which used type B ol 1 , experienced an average life of 300 hours, and subsur- 
face fatigue failures In the form of pitting. The type B oil was a synthetic 
paraffinic oil with five volume percent of an extreme-pressure (EP) additive. 
This EP additive contains sulfur and phosphorus to enable the oil to keep an 
adequate EHD film thickness, minimizing metal to metal contact, even under 
extreme pressure. 

The failure modes experienced by the various gears can be categorized Into 
four basic damage related groups. These groups are: 1) Heavy wear and scor- 

ing, 2) Tooth breakage, 3) Single pits, and 4) Distributed pitting. The groups 
are based on the types and magnitude of damage found on the gears at the end of 
their runs. As Indicated In Table I, run 8 was not classified as failed 
because It ran the maximum of 520 hours without exceeding the rig's vibration 
limit. The point of failure on the fatigue rig was determined by an overall 
vibration level from an accelerometer mounted on the rig. When this level 
reached the preset threshold, the rig shutdown and the run was classified as 

failed. Because the gears of run 8 experienced damage similar to another run. 

It Is Included In that run's group. A more detailed description Is given In 
the following sections. 

HEAVY WEAR AND SCORING: 

Runs 1, 2, and 4 through 7 all experienced heavy wear and scoring uni- 
formly over all the gear teeth. Figures 3.2 and 3.3 Illustrate the wear pat- 
terns observed. As seen in these figures, the driven gear teeth experienced 
heavy wear and scoring at the tip portion of the tooth, a band of light pitting 

at the double to single tooth contact point, and heavy wear In the root region 

of the tooth. The driver gear teeth experienced similar wear patterns to the 
driven gear teeth, with the exception that the top portion of the tooth exper- 
ienced lighter wear and scoring. Run 3 experienced the same wear patterns, 
although to a lesser degree due to the premature tooth failure at 9 hours Into 
the test. The tooth surfaces appeared extremely course, especially In the 
regions of heavy wear and scoring. 

TOOTH BREAKAGE 

Runs 3 and 7 experienced tooth breakage during the fatigue test conducted. 
One tooth broke during run 3, causing the rig to shut down. A photograph of 
the fracture surface on the driven gear In run 3 Is shown In figure 3.4. No 
detailed tooth crack propagation analysis was done on this tooth, however by 
visual inspection, there appeared to be only three distinct bands in the frac- 
ture area, with little evidence of fracture surface smoothing. This indicates 
that the time between crack formation and tooth fracture was relatively short. 
Run 7 also experienced tooth fracture; however, a failure in the automatic shut 
down mechanism did not stop the rig soon enough to prevent extensive damage to 
both gears. Thus no tooth fracture observations could be performed on run 7 
gears . 
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SINGLE PITS 

Runs 8 and 9 experienced similar damage In the form of single pits. 

Flqures 3.5 and 3.6 Illustrate the wear pattern found on the gears In run 8 at 
the end of the run. As seen In these figures all but one of the teeth on the 
driver and driven gear appear to have no sign of wear or pitting. A moderate 
pit appears on one tooth on the driver gear, and a larger pit appears on one 
tooth of the driven gear. Figures 3.7 and 3.8 Illustrate the condition of the 
qears in run 9 at the point of failure. As seen In these figures, run 9 
appears to have a similar failure pattern as found in run 8, with one excep- 
tion. The large pit on the driven gear In run 9, although roughly 
size as the large pit In run 8, has It's major axis orientated across the gear 
tooth rather than along the tooth surface, as found In run 8. The 
found on the gears In runs 8 and 9 are characterized as single tooth defects. 

Runs 10 and 11 experienced similar damage In the form of distributed pit- 
ting. Figures 3.9 and 3.10 Illustrate the wear pattern present on the gear 
teeth in run 10. As seen in these figures, the driver gear exhibits * 
pattern of pitting bands with various Intensities over 60 percent of the teeth. 
The driven gear displayed a similar pattern over a smaller region of teeth. 
Figures 3.11 and 3.12 Illustrate the wear pattern Present on the gear teeth of 
the gears used In run 11. As seen In these figures, the distributed pitting 
bands appear In a similar fashion to those In run 10, although to a lesser 
extent. 


3.2 INSTRUMENTATION SET-UP 

Flqure 3 13 illustrates the test Instrumentation set-up used In this 
study. As seen In this figure, two accelerometers. PCB Model 303A, were 
mounted on the rig close to the test gear mesh. The vibration signal from 
these accelerometers were recorded on a high precision tape recorder, g 
Heston Sabre VII, along with a once per revolution signal and a time code si g 
nal . A one minute signature was collected every three hours by using a timer 
to control recording time. The once per revolution signal was provided by a 
photon sensor that produced a narrow time pulse (.202 msec) for each revolution 
of the shaft. This signal was used for time synchronous averaging, and for 
determining the actual rotational speed. The time code signal was used during 
tape playback to distinguish between the separate data points, and to provide 
the exact time and day for each data point recorded. ... , ,. h . 

After the runs were recorded, the data was then analyzed by replaying the 
tape into a single channel dynamic signal analyzer, with a dynamic range of 
80 dB, and transferring it to a personal computer. The raw data was digitized 
and averaged by the analyzer and sent via a general Purpose Interface bus to a 
IBM compatible personal computer. The PC used ;. s *JJ da J d ?EK-488 interface 
board to control the analyzer and download the digitized data to DOS f lies on 
the PC's disc drive. The data transferred was the averaged amplitude and phase 
portions of the Fourier transform. The data was then analyzed by using several 
FORTRAN programs, each based on the various prediction techniques inve ^^ated. 

The primary mesh frequency was an important factor In choosing and s 9 
the test instrumentation. For the 28 teeth test gears operating at 10,000 rpm, 
the primary and first harmonic of the mesh frequency occurs at a ppr°x 1 ma t e 1 y 
4700 Hz and 9400 Hz, respectively. The accelerometers were chosen wi th a fre- 
quency range of 0 to 10,000 Hz to capture the mesh frequencies. The recorder- 
tape speedwas set at 762 mm/sec (30 in/sec) with wideband group I FM recording 
electronics to give a recording frequency range of 0 to 20 000 Hz, with a 
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dynamic range of 50 dB. FM recording was used to give good amplitude sta- 
bility, Independent of the storage characteristics of the recording tape. 

Using a tape reel with a 2800 m length (9200 ft), one reel translated to 
184 hours of coverage, or 7.7 days. This was considered a reasonable time 
interval, since some of the tests ran continuously for nearly 21 days. The 
analyzer was set up to time average the signal over 40 averages for each data 
point, using the shaft synchronous signal as the input trigger. The FFT was 
performed over the range 0 to 12,500 Hz. 

During the first stages of the data analysis process, It was discovered 
that the side accelerometer provided the same Information as the top accelerom- 
eter. In fact, the meshing frequency components In the side accelerometer sig- 
nal were not as dominant as In the top accelerometer signal. This is probably 
due to the fact that the top accelerometer Is orientated to sense the vertical 
vibrations of the box. This is coincident with the direction of the primary 
component of the tooth mesh loads, resulting from the side by side arrangement 
of the test gears. Thus, the vibration data from the top accelerometer was 
used primarily In this analysis, with the side accelerometer data used as a 
backup. 

3.3 FREQUENCY RESPONSE MEASUREMENTS OF FATIGUE RIG 

During the Initial stages of the data analysis, It was discovered that 
the recorded time signals contained frequency data not coherent with the gear 
meshing frequencies. Figures 3.14 through 3.17 Illustrate the presence of the 

extraneous frequencies as run 7 progressed from Initiation to failure. This 

phenomenon was found in all of the runs recorded. As seen In these figures, 
certain frequencies, not coincident with the gear mesh frequencies or asso- 
ciated sidebands, Increase drastically through the run cycle. The primary mesh 
frequency, highlighted by an arrow In the figures. Is In most cases lower In 
magnitude than the extraneous frequencies. The amplitudes near 3000 and 

6000 Hz appear to be the most dominant. If the gear failure prediction tech- 

niques were applied to the signals given In Figures 3.14 through 3.17, results 
not Indicative of the actual gear condition would dominate. Therefore, It was 
necessary to Investigate these extraneous frequency components and determine 
their origin. If these components could be attributed to a known source, they 
then could be removed from the signal with confidence. 

A number of possible sources of the extraneous frequency components were 
Investigated. Because the signal is synchronous time averaged with the gear 
shaft rotation, the only possible periodic sources coincident with the shaft 
speed are the roller bearings, ball bearings, and the slave gear. The location 
of these components are given In Figure 3.1. Based on the assumptions of a 
single defect and that the rolling elements remain in rolling contact, bearing 
defect frequencies were calculated using existing equations. [11] These fre- 
quencies covered outer race, inner race, and rolling element defects, and fell 
within the 64 to 1400 Hz range. These explain some lower frequency components, 
but they were too low for the dominant frequencies In question. The slave gear 
primary mesh frequency was found to be 6084 Hz, which Is near some of the 
extraneous frequencies, but it was not exactly coincident. The next step in 
Identifying the frequencies In question was to determine the modal parameters 
of the rig using frequency response measurements. 

Frequency response measurements performed on the fatigue rig Illustrated 
how the dynamics of a transfer path can critically alter a signal. The fre- 
quency response measurements were taken at two locations on the box with the 
reference point being the top accelerometer location, as illustrated in 
Figure 3.18. For these tests, the impact method was used. A modal hammer with 
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a steel tip and force transducer was used as the Input to the system, with the 
top accelerometer serving as the output signal. The force transducer and 
accelerometer were Input Into a two channel dynamic signal analyzer capable of 
calculating the frequency response function. Twenty averages were taken for 
each measurement. The FFT of the Input Impulse Is given In Figure 3.19. As 
seen In this figure, the Impulse adequately excites the frequencies in 
question. 

Results from the frequency response measurements directly correlated the 
unexplained frequency components to the path modes. The frequency response 
measurement of the path from the top of the bearing cover to the top accelerom- 
eter location is given in Figure 3.20. The coherence function for this meas- 
urement, shown In Figure 3.21, Indicates the output Is almost totally a func- 
tion of the Input up to 10,000 Hz. Figures 3.22 and 3.23 give the frequency 
response measurement and Its' coherence for the path from the gear shaft to the 
top accelerometer location. As seen in Figures 3.20 and 3.22, the vibration 
path from the gear to the top accelerometer drastically alters the signal due 
to the natural frequencies Inherent in the box. As the gears experience wear, 
the meshing energy redistributes from the primary meshing frequencies to a wide 
range of frequencies in the frequency spectrum. When these frequencies become 
coincident with the natural frequencies associated with the transfer path the 
corresponding amplitudes Increase. This phenomenon Is clearly Illustrated in 
the 3000 Hz region and in the region within the 5000 to 8000 Hz band in 
Figures 3.14 through 3.17. The frequencies In these regions are coincident 
with the natural frequencies associated with the transfer path, as evident In 
Figures 3.22 and 3.23. 

Most of the frequencies in question were found to be a direct result of 
the structural dynamics of the transfer path and as such are not desired as 
inputs for the various failure prediction techniques. Therefore these frequen- 
cies were removed from the processed signals by filtering techniques. The 
functions illustrated in Figures 3.20 and 3.22 were used to determine the fre- 
quencies to remove. Careful attention was placed on eliminating only those 
frequencies known to be dominant, and to allow frequencies adjacent to the mesh 
frequencies to be passed. Fortunately, the primary mesh frequency region had 
no dominant natural frequencies present; however some of the upper sidebands 
were eliminated with the filtering. All of the analysis methods were applied 
to the filtered signal. It is acknowledged that some of the desired signal was 
eliminated with the filtering; however, It is assumed that the errors associated 
with the filtering process are less than those associated with leaving the 
dynamic effects of the transfer path In the signal. The frequency bands most 
commonly filtered out were between 400 and 1500 Hz, 2400 and 3700 Hz, and 5200 
and 8400 Hz. 
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Figure 3.1 : Gear Fatigue Test Apparatus. 

a) Schematic Diagram. 

b) Cutaway View. 
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Figure 3.2 : Illustration of Tooth Damage Found on Gears in Runs 1, 

2, and 4 through 6. 
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b) driven gear teeth 

Figure 3.3 : Photographs of uniform damage found 
on gears in Run 1 
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Figure 


.4 : Photograph of tooth fracture area on 
driven gear In Run 3 
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Figure 3.5 : Illustration of Tooth Damage Found on Gears in Run 8. 
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a) driver gear tooth 



b) driven gear tooth 

Figure 3,6 : Photographs of tooth damage found 

on gears in Run 8 
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Figure 3.7 : Illustration of Tooth Damage Found on Gears in Run 9. 
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b) driven gear tooth 

Figure 3.8 : Photographs of tooth damage found 
on gears in Run 9 
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Figure 3.9 : Illustration of Tooth Damage Found on Gears in Run 10. 
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a) Tooth # 12, driver gear 



b) Tooth # 8, driver gear 
Figure 3*10 : Photographs of tooth damage found 

on gears in Run 10 
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d) Tooth # 10, driven gear 

Figure 3.10 : Concluded. 
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Figure 3.11: Illustration of Tooth Damage Found on Gears in Run 11. 
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b) Tooth # 16, driver gear 
Figure 3.12 : Photographs of tooth damage found 
on gears in Run 11 
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d) Tooth # 4, driven gear 

Figure 3.12 : Concluded. 
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Figure 3.13: Illustration of Test Set-up. 










Figure 3,14: Frequency Spectrum of Run 7 at Start of Run. 



Figure 3.15: Frequency Spectrum of Run 7 After 32 Hours of Running. 
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Frequency Spectrum of Run 7 After 44 Hours of Running. 
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Figure 3.17 : Frequency Spectrum of Run 7 After 50 Hours of Running. 
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Figure 3.18 : Illustration of Frequency Response Measurement 

Locations . 
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Figure 3.19 : Power Spectrum of Input to Frequency Response 

Measurements . 
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Figure 3.20 : Magnitude of Frequency Response Function, Test 1, 

Between the Top Accelerometer Location and the Top of 
the Bearing Housing. 



Frequency (Hz) 

Figure 3.21 : Coherence Function of Test 1 Measurements, Between the 

Top Accelerometer Location and the Top of the Bearing 
Housing. 
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CHAPTER IV 

APPLICATION AND RESULTS 


The main portion of this study was the application of the various gear 
fault methods to experimental data. Chapter 4 briefly explains the various 
computer programs developed that apply each of the techniques. Also, this 
chapter presents and discusses the results of these techniques. 

4.0 PROGRAMS DEVELOPED 

To apply the various predictive techniques to the experimental data, sev- 
eral computer programs were developed. The programs are written In FORTRAN, 
and listed In the appendices. For those methods requiring Fourier analysis, 
the standard FFT algorithm developed by Cooley and Tukey [12] Is Included in 
each routine. The computer programs use the data files created using the 
dynamic signal analyzer. The results of the programs are stored In data files 
that can be plotted using commercially available routines. 

The programs use the various equations and theories presented In Chapter 
II. Program EFMO. FOR, listed In Appendix A, calculates the FMO parameter for 
each time data point In the run. Program TFM4.FOR, listed In Appendix B, 
applies the FM4 techniques to the data by calculating the normalized kurtosls 
and standard deviation values of the difference file for each time data point. 
Program HILB. FOR, listed In Appendix C, uses the Hilbert transform technique 
to calculate the Instantaneous phase variations as a function of shaft posi- 
tion for each time data point. Program PARAM.FOR, listed In Appendix D, calcu- 
lates the crest factor, sideband level factor, and energy ratio for each time 
data point. 

The computer programs were verified by Inputting known functions. As an 
example of this, the kurtosls and standard deviation routines In program 
TFM4.F0R were verified by inputting a sine wave and square wave, and comparing 
the results with exact solutions. Table II illustrates the results of this 
verification test. The square wave did not provide as high a correlation as the 
sine wave, mainly because the input square wave contained some irregularities. 


TABLE II: Program TFM4. FOR routines verification 


Routine - Input 

Measured value 

Actual Value 

1 Difference 

Standard Deviation 

- Sine wave 

- Square Wave 

1.3472 

1.9180 

1.3435 
1 . 9000 

0.28 % 
0.95 % 

Normalized Kurtosls 

- Sine Wave 

- Square Wave 

1.4896 
1 .0404 

1 . 5000 
1 .0000 

0.69 t 
4.04 X 
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4.1 RESULTS 


FMO METHOD , ,, 

The FMO method did very well at detecting a majority of the heavy wear 
and scoring damage experienced by runs 1, 2, and 4 through 6, as It was 
designed to do. Figures 4.1 through 4.5 plot the three parameters: a) FMO, 

b) Peak-to-peak level (the FMO numerator), and c) Sum of mesh amplitudes (the 
FMO denominator), for these runs. As seen in these figures, all of the runs, 
except run 1, exhibit an increase In the FMO parameter, plot a, with Increasing 
run time. These same runs also show a decrease In the primary meshing frequency 
and second harmonic amplitude sum, plot c, with increasing run time. It Is 
Interesting to note that in Figure 4.3, run 4, the FMO plot Indicates heavy 
damage at 49 hours into the test, due to a decrease In meshing frequency ampli- 
tudes, whereas the peak-to-peak level, plot b, starts Increasing at 55 hours 
and peaks at 64 hours. Run 1 did not provide good FMO trends even though the 
gear experienced similar damage. 

FMO is designed to detect major tooth faults such as the breakage exper- 
ienced during runs 3 and 7. Because no data was collected at the time of break- 
age, or immediately after it occurred, FMO was unable to be applied to these 
runs . 

FMO was applied to runs 8 and 9, which experienced single tooth pits, even 
though FMO is not designed to detect single tooth faults. As seen in Figure 
4.6, run 8 shows a gradual overall increase in FMO values with run time, and a 
gradual overall decrease In meshing frequency amplitudes. This trend In run 8 
could be attributed to its long run time of 520 hours, which may have resulted 
In uniform wear not observable by visual Inspection. Run 9 showed no logical 

trend with FMO. , ... 

FMO was also applied to runs 10 and 11, where the gears experienced dis- 
tributed pitting damage. As seen in Figures 4.7 and 4.8, the FMO values for 
both runs Increase near the end of the run time. The meshing frequency ampli- 
tudes show an overall decrease with increasing run time, Indicating the exist- 
ence of heavy wear. Since FMO does not respond to specific tooth damage, it 
is reasonable to assume that the pitting occurred over enough teeth to act as 
a uniform wear phenomenon, and thus was capable of being detected by FMO. It 
Is Interesting to note that for run 10, Figure 4.7, the FMO value Increased 
sharply at the 220 hours point. At this same point In time, the meshing fre- 
quencies amplitude sum decreased sharply, possibly Indicating the time when 
most of the major distributed pitting happened. The peak-to-peak level at 
this period contained no major fluctuation, only a gradual Increase in level. 

FM4 METHOD 

The normalized kurtosis parameter was unable to detect the single tooth 
faults found in the test gears. Runs 8 and 9 experienced specific faults in 
the form of one large and one small pit. As seen in Figures 4.9 and 4.10, 
these defects had limited effect on the kurtosis parameter. The fatigue cracks 
that resulted in broken teeth in runs 3 and 7 were also undetected by the 
kurtosis parameter, as seen in Figure 4.11 for run 7. The normalized kurtosis 
values for the two data points of run 3, not plotted, were 2.98 and 3.04. It 
should be noted that the data points collected at the end of runs 3 and 7 were 
taken 3 and 2 1/2 hours, respectively, prior to the point of tooth fracture. 

Only run 11 had limited success using the normalized kurtosis parameter. As 
seen in Figure 4.12, a normalized kurtosis value of 5.3 was registered at 
88 hours of operation; however, the value decreased to approximately 3.0 at 
96 hours. One explanation for this trend is that only one large pit was 
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present at the 88 hour mark. As more pits formed, the number of peaks 
Increased causing the normalized kurtosls value to decrease. 

The standard deviation parameter, although used In FM4 for single tooth 
failure detection, proved to be a good Indicator of heavy wear. Runs 1, 2 and 
4 through 6 all experienced heavy wear and scoring. The standard deviation 
plots for these runs, given In Figures 4.13 through 4.17, all show clear trends 
of Increasing values near the end of each run with only minor fluctuations. 

The standard deviation plot of run 1, Figure 4.13, Indicates a definite trend 
as compared to the FMO plot of run 1, which did not respond to the heavy wear. 

Plots of the squared difference signal of FM4 provided some general Infor- 
mation on the wear condition of the tooth surfaces. The plots of the squared 
difference signal are presented as a function of gear rotation. Two plots are 
given In each figure. The top graph (a) represents the data collected at the 
beginning of the run and the bottom graph (b) represents the data collected at 
the end of the run. The runs that experienced heavy wear and scoring show def- 
inite Increases In the difference signal between the start and end plots 
throughout the gear rotation. An example of this Is given In Figure 4.18. The 
specific pits of runs 8 and 9 were not seen in their squared difference signal 
plots. The wear pattern on the driver gear In run 10 Is reflected In the 
squared difference signal given In Figure 4.19. Run 11 experienced similar 
wear as run 10; however, the wear pattern was not easily detected in Its 
squared difference signal plot. 

HILBERT TRANSFORM METHOD 

The Hilbert transform method was primarily developed to detect fatigue 
cracks using the phase modulation function described In Section 2.1. Runs 3 
and 7 are the only runs that experienced tooth fracture due to probable fatigue 
cracks. Figure 4.20 plots the phase modulation function for the two data time 
Intervals of run 3. Several phase shifts can be seen In the first plot, 2 
hours Into the run; however, they are not reflected In the second plot repre- 
senting the last data point. Nothing in the last data point plot Indicates 
the presence of a fatigue crack, i.e. no large phase lags present. Figure 4.21 
plots the phase modulation function for the last two data time Intervals of 
run 7. The only possible Indication of a fatigue crack starting Is the phase 
lag at approximately half way Into the shaft rotation, near the 180 degrees 
point. The phase lag starts during the second to last time interval, 47 hours 
into the run, and grows to the size seen In the last Interval, 50 hours Into 
the run. More data after this last data point is required to claim with any 
certainty that this phase lag represents an actual fatigue crack. Again it 
must be noted that the last data points for runs 3 and 7 were taken 3 and 2 
1/2 hours, respectively, before tooth fracture. 

CREST FACTOR AND SIDEBAND LEVEL FACTOR 

Both the crest factor and sideband level factor are designed to respond 
to signals with Impulsive vibration sources, specifically tooth breakage. 

Runs 3 and 7 were the only runs that experienced tooth fracture. The crest 
factor and sideband level factor for run 7 are plotted In Figure 4.22. Run 3 
was not plotted, since it had only two time Intervals recorded before fracture 
occurred. As seen In Figure 4.22, neither the crest factor (plot a) or the 
sideband level factor (plot b) display any indication that a tooth fracture 
was going to occur. These parameters may be sensitive to tooth breakage only 
after it has happened. Unfortunately, no data was collected during or after 
the fracture occurred. 
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ENERGY RATIO ,, 

The energy ratio is designed to be a robust Indicator of heavy wear. 

Runs 1, 2, and 4 through 6 experienced heavy wear and scoring. The energy 
ratio graphs of these runs are given in Figures 4.23 through 4.27. As seen in 
these figures, the energy ratio does not provide as good an indication of wear 
as the FMO and the standard deviation methods. The energy ratio graphs for 
runs 10 and 11 are given In Figure 4.28 and 4.29, respectively. Of these two 
runs, only run 10 had an Increase in its' energy ratio parameter. 
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Figure 4.1 : Plot of FMO Parameters for Run 1 

a) FMO Parameter. 

b) Numerator of FMO. 

c) Denominator of FMO. 
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Figure 4.2 : Plot of FMO Parameters for Run 2 

a) FMO Parameter. 

b) Numerator of FMO. 

c) Denominator of FMO. 
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Figure 4.5 : Plot of FMO Parameters for Run 6 

a) FMO Parameter. 

b) Numerator of FMO. 

c) Denominator of FMO. 
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Figure 4.6 : Plot of FMO Parameters for Run 8. 

a) FMO Parameter. 

b) Numerator of FMO. 

c) Denominator of FMO. 
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Figure 4.8 : Plot of FMO Parameters for Run 11. 

a) FMO Parameter. 

b) Numerator of FMO. 

c) Denominator of FMO. 
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Figure 4.9 : Plot of Normalized Kurtosis (FM4) for Run 8 
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Figure 4.10 : Plot of Normalized Kurtosis (FM4) for Run 9 
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Figure 4.11 : Plot of Normalized Kurtosis (FM4) for Run 7. 



Figure 4.12 : Plot of Normalized Kurtosis (FM4) for Run 11. 
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Figure 4.13 : Plot of Standard Deviation of Difference Signal (FM4) 

for Run 1. 



Figure 4.14 : Plot of Standard Deviation of Difference Signal (FM4) 

for Run 2. 
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jure 4.15 : Plot of Standard Deviation of Difference Signal (FM4) 

for Run 4. 
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Figure 4.16 : Plot of Standard Deviation of Difference Signal (FM4) 

for Run 5. 
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Figure 4.17 : Plot of Standard Deviation of Difference Signal (FM4) 

for Run 6. 
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Figure 4.18 : Plot of Squared Difference Signal For Run 4. 

a) At Start of Run. 

b) At End of Run. 
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Figure 4.19 : Plot of Squared Difference Signal For Run 10. 

a) At Start of Run. 

b) At End of Run. 
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Figure 4.20 : Plot of Phase Modulation Function for Run 3. 

a) At a Run Time of 2 Hours. 

b) At a Run Time of 5 Hours. 
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Figure 4.21 : Plot of Phase Modulation Function for Run 7 

a) At a Run Time of 47 Hours. 

b) At a Run Time of 50 Hours. 
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Figure 4.22 : Plot of Crest Factor and Sideband Level for Run 7 

a) Crest Factor. 

b) Sideband Level Factor. 
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Figure 4.23 : Energy Ratio Plot for Run 1 
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Figure 4.24 : Energy Ratio Plot for Run 2 
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Figure 4.25 : Energy Ratio Plot for Run 4 
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Figure 4.26 : Energy Ratio Plot for Run 5 
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Figure 4.27 : Energy Ratio Plot for Run 6. 
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Figure 4.28 : Energy Ratio Plot for Run 10. 
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Figure 4.29 : Energy Ratio Plot for Run 11 



CHAPTER V 

CONCLUSIONS AND RECOMMENDATIONS 


The type and extent of the damage found on the eleven test gears can be 
classified Into four major failure modes. The first mode can be described as 
heavy wear and scoring (runs 1,2 and 4 through 6). The second mode Is tooth 
breakage (runs 3 and 7). The third failure mode is a result of single pits 
(runs 8 and 9). And the last mode is described as distributed pitting (runs 
10 and 11). These classifications will be used to evaluate the overall per- 
formance of the various methods. 

HEAVY WEAR AND SCORING 

The FMO parameter and the standard deviation of the difference signal 
from FM4 did very well at detecting this heavy wear condition. Results from 
both techniques support the theory that as a gear wears the meshing energy 
redistributes from the meshing frequencies to Its sidebands and beyond. 

The energy ratio method does not predict wear as well as the FMO and the 
standard deviation methods. The most probable reason for this is in the denom- 
inator of the energy ratio. The denominator is the standard deviation of the 
regular signal, or the meshing frequencies and their first order sidebands. 

The first order sidebands should not be considered part of the regular signal 
for this parameter. This is because they Increase In amplitude similar to the 
higher order sidebands as the gear wears. 

An enhanced method using a combination of these techniques could prove to 
be a reliable uniform gear wear detection technique. One such method could 
use the standard deviation of the difference signal divided by the sum of the 
amplitudes of the meshing frequencies (primary and first harmonic). Because 
wear becomes detectable in the overall vibration level only when It reaches a 
critical stage, a robust wear Indication parameter would be highly useful. 

TOOTH BREAKAGE 

The fatigue cracks that resulted in the broken teeth in runs 3 and 7 were 
not detected by any of the methods. The normalized kurtosis parameter of FM4, 
the phase modulation function from the Hilbert transform technique, the crest 
factor, and the sideband level factor are all conditioned to react to tooth 
cracks to varying degrees; however, none detected any fault prior to tooth 
fracture. The most probable reason for this can be due to the high speed and 
high loading conditions the test gears are subjected to. With these operating 
conditions, the time elapsed between initiation of the fatigue crack and even- 
tual tooth fracture is probably orders of magnitude lower than the three hour 
data interval time. Since the last time records collected from runs 3 and 7 
were 3 and 2 1/2 hours before actual tooth fracture, respectively, data 
necessary to indicate the fatigue cracks were missed. 

SINGLE PITS (TOOTH SPECIFIC FAULTS) 

The single pits were not detected by the normalized kurtosis in the FM4 
technique. The normalized kurtosis of the difference signal is designed to 
detect tooth specific faults; however, no indication of these defects were 
recorded with this parameter. Although they physically appear large, they may 
have not affected the signal to the point detectable by FM4. 

DISTRIBUTED PITTING 

The distributed pitting damage was detected with the the FMO parameter. 

It Is theorized that the pitting happened over enough of the teeth to act as a 
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uniform wear phenomenon, thus becoming detectable to FMO. The square of the 
difference signal from the FM4 method did appear to reflect the actual wear 
pattern on some of the runs. 

RECOMMENDATIONS 

It Is recommended that any future tests of this nature be Implemented on 
an on-line basis for three major reasons. The first stems from the fact that 
a passive system cannot call attention to events which happen faster than the 
time between the data time intervals. A continuous, or on-line, system can 
monitor the test and record data as fast as the electronics will allow when 
certain parameters signal a significant event in progress. This may be the 
only way to provide data on the fatigue cracks that lead to tooth fracture on 
this rig. Another major reason is that in an on-line system it is possible to 
have the system stop the rig when a parameter detects a fault. The fault can 
be analyzed to determine its extent and nature at the exact time it is 
detected. This would provide an excellent means of verifying and linking cer- 
tain parameters with certain faults as they naturally happen. The final advan- 
tage of an on-line system Is the ability of storing data in digital form, 
eliminating the FM recorder. In this test, because the dynamic range of the 
recorder was only 50 dB, some of the lower energy components of the signal were 
lost. Using the analyzer to digitize and store the data directly would 
Increase the dynamic range to 80 dB, thereby capturing more of the lower energy 
signal s . 

It is also recommended that frequency response measurements be performed 
between the gear shaft and the proposed measurement sites. The dynamics of the 
transfer path In this study were found to significantly affect the output sig- 
nal, especially as the gear began to wear. An understanding of the structural 
dynamics associated with the transfer path is required to apply prediction 
techniques only to the portion of the signal reflecting the actual gear mesh 
dynamics. 
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APPENDIX A - SOURCE CODE FOR EFMO. FOR 


c c 
c c 

C PROGRAM: EFMO . FOR CREATED: 05-05-89 REVISED: XX-XX-XXXX C 
C C 
C PROGRAMMER: J.J. ZAKRAJSEK C 
C C 
C THIS PROGRAM COMPUTES AN ENHANCED VERSION OF THE COURSE FAULT C 
C DETECTION PARAMETER (EFMO). IT IS THE BASIC FMO PARAMETER, C 
C AS PROPOSED BY STEWART IN HIS 1977 PAPER, HOWEVER, THE PEAK C 
C TO PEAK VALUE USED IN EFMO IS THE PEAK TO PEAK OF THE C 
C RECONSTRUCTED TIME SIGNAL AFTER BAND WIDTH FILTERING IS C 
C PERFORMED. THE RECONSTRUCTED TIME SIGNAL USED IS THE "FUME" C 
C FILE CREATED AFTER A-B-C FILTERING DURING FM4 CALCULATIONS IN C 
C PROGRAM TFM4 . FOR . THE AMPLITUDES AT THE MESH FREQUENCY AND ITS C 
C HARMONIC ARE RETRIEVED FROM PREVIOUS FMO CALCULATIONS. C 
C C 
C C 


c 

DIMENSION TIME(50) , ECFDP(50) , PKTPK(50) ,ASUM(50) ,FSIG(220) 
CHARACTER*22 FMOF , EFMOF , LIST2 , FTIMEF 
C 

C ENTER INITIAL DATA AND FILE NAMES 
C 

PRINT *, 'ENTER NUMBER OF DATA POINTS' 

READ(*,*) K 

PRINT *, 'ENTER OUTPUT DATA FILE NAME (\TXXE*FM0 . DAT) ' 

PRINT * , ' ( * - FILTER CONFIG USED IN RECONSTRUCTED FILTER) ' 

READ (*, ' (A22) ' ) EFMOF 

OPEN (7 , FILE=EFMOF, IOSTAT=IOERR) 

IF (IOERR .GT. 0) GOTO 999 

PRINT *, 'ENTER FILE NAMES ARRAY FILE (\TXXLIST2 . FIL) ' 

READ(* , ' (A22) ' ) LIST2 

OPEN (8 , FILE-LIST2 , IOSTAT=IOERR) 

IF(IOERR .GT. 0) GOTO 999 

PRINT *, 'ENTER PREVIOUS TFMO FILE (\TXXFMO.DAT)' 

READ(*, ' (A22) ') FMOF 

OPEN (9 , FILE=FM0F, IOSTAT=IOERR) 

IF( IOERR . GT. 0) GOTO 999 

PRINT *, 'ENTER NUMBER OF TIME DATA POINTS PER CYCLE (INTEGER)' 
READ(* , *) NCYC 
C 

C START TIME LOOPING, PROCESSING DATA AT EACH TIME INTERVAL 
C 

DO 300, J-1,K 

PRINT *, 'ACCESSING AND PROCESSING DATA FILE NO:',J 
READ(8 , 800) FTIMEF 
800 FORMAT (A18) 

READ (9, 900) TIME(J) ,ASUM(J) 

900 FORMAT(F8 . 2 , 44X , F12 . 8) 

OPEN (4 , FILE-FTIMEF , IOSTAT=IOERR) 

IF (IOERR .GT. 0) GOTO 999 
C 

C READ IN AND FIND THE MAX AND MIN VALUES OF THE FILTERED TIME RECORD 
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C 

AMAX - 0.0 
BMAX =0.0 
DO 100, I-l.NCYC 
READ (4 , 400) FSIG(I) 

400 FORMAT (20X.F16. 8) 

IF( FSIG(I) .GT, AMAX ) THEN 
AMAX = FSIG(I) 

END IF 

IF( FSIG(I) .LT. BMAX ) THEN 
BMAX = FSIG(I) 

END IF 

100 CONTINUE 

C 

C DETERMINE PK TO PK OF FILTERED SIGNAL, AND ENHANCED FMO VALUE 
C 

PKTPK(J) - AMAX - BMAX 
ECFDP(J) = PKTPK ( J ) /ASUM ( J ) 

C 

C CLOSE LOOP DEPENDANT FILE, AND WRITE CALAULTED DATA 
C 

CLOSE (4 , IOSTAT=IOERR) 

IF (IOERR .GT. 0) GOTO 999 

WRITE (7,700) TIME(J) ,ECFDP(J) ,PKTPK(J) 

700 FORMAT(F8 . 2 , F10 . 4 , F10 .4) 

300 CONTINUE 
C 

PRINT * , ' EFMO VALUES STORED IN FILE ; ' , EFMOF 
STOP 

999 PRINT*, 'ERROR NO. ', IOERR 
END 
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APPENDIX B - SOURCE CODE FOR TFM4. FOR 


c c 
c c 

C PROGRAM: TFM4 . FOR REVISION: NEW C 
C REVISION DATE: NEW C 
C PROGRAMMER: J.J. Zakrajsek DATE CREATED: 04-13-89 C 
C C 
C THIS PROGRAM COMPUTES THE FAULT DETECTION PARAMETERS AND THE C 
C RELATED MATHEMATICS FOR THE FM4 DIAGNOSTIC THEORY PROPOSED C 
C IN STEWART'S 1977 PAPER " SOME USEFULL DATA ANALYSIS C 
C TECHNIQUES FOR GEARBOX DIAGNOSTICS." THESE INCLUDE PARAMETERS C 
C SUCH AS STANDARD DEVIATION AND KURTOSIS. C 
C C 
C C 


£****■***************************■*************************•***■******■***£ 

C 

DIMENSION TIME(50) , DYO(220) , DSYO(220) ,RMAG(401) 

DIMENSION FREQ(401) , FMAGRMS(401) , FMAG(401) ,PHAS(401) ,FFMAG(401) 
DIMENSION TO(1024) ,YO(1024) ,FYO(1024) ,RFYO(1024) 

CHARACT ER* 2 2 TIMEF , FM4F , LIST1 , LIST2 , AMPLF , PHASEF , ATIME 
CHARACTER* 18 FTIME , RTIME , NTIME , STIME 
COMPLEX X(1024) 

C 

C ENTER GENERAL PARAMETERS AND FILE NAMES 
C 

PRINT *, 'ENTER NUMBER OF DATA POINTS (INTEGER)' 

READ(*,*) K 

PRINT *, 'ENTER MESHING FREQUENCY, (XXXX.XX) ' 

READ ( * , * ) PRIM 

PRINT *, 'ENTER ROTATIONAL SPEED IN REV/SEC (Hz)' 

READ(*,*) SPD 

PRINT *, 'ENTER NUMBER OF TIME DATA POINTS PER CYCLE' 

PRINT *, '( INTEGER, 220 MAXIMUM )' 

READ (*,*) NCYC 
C 

C DEFINITION OF OPTIONAL FILTER BANDS A, B, AND C 
C 

PRINT *, '3 BAND WIDTH DIGITAL FILTERS CAN BE SELECTED.' 

PRINT *, 'DATA EXISTING WITHIN THE THREE BANDS WILL BE FILTERED' 
PRINT *, 'OUT OF THE ANALYSIS BEFORE IT BEGINS FM4 CALCULATIONS' 
PRINT *, ' 

PRINT *, 'TO USE FILTERS ENTER (1)' 

PRINT *, 'TO SKIP FILTERS ENTER (0)' 

C 

READ (*,*) NN 
IF (NN .LT. 1) GOTO 123 
C 

PRINT *, 'ENTER STARTING FREQUENCY OF FILTER A' 

READ ( * , * ) AFILTS 

PRINT *, '*************** ENTER END FREQUENCY OF FILTER A' 

READ(* , *) AFILTE 

PRINT *, 'ENTER STARTING FREQUENCY OF FILTER B' 

READ(* , *) BFILTS 

PRINT *, '*************** ENTER END FREQUENCY OF FILTER B' 

READ ( * , * ) BFILTE 
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PRINT *, 'ENTER STARTING FREQUENCY OF FILTER C' 

READ(*,*) CFILTS 

PRINT *, **************** ENTER END FREQUENCY OF FILTER C' 
READ(*,*) CFILTE 
123 CONTINUE 


C 

C 

C 


c 

C 

c 


INPUT PRIMARY FILE NAMES FILES 


PRINT *, 'ENTER 
READ(* , '(A22) ') 
PRINT * , ' ENTER 
PRINT * , ' ( * - 
READ(* , ' (A22) ') 
PRINT *, 'ENTER 
READ(* , ' (A22) ') 
PRINT * , ' ENTER 
READ(*, ' (A22) ' ) 


TXXTIME . FIL 
TIMEF 

TXXFM4* . FIL 
FILTERS A-B-C 
FM4F 

TXXLIST1 . FIL 
LIST1 

TXXLIST2 , FIL 
LIST2 


INPUT FILE NAME' 

OUTPUT FILE NAME' 
VERSION) ' 

INPUT FILES FILENAME' 
OUTPUT FILES FILENAME' 


OPEN PRIMARY FILES 


OPEN ( 1 , FILE=TIMEF , IOSTAT-IOERR) 
IF(IOERR .GT. 0) GOTO 999 
OPEN ( 2 , FILE-FM4F , IOSTAT-IOERR) 
IF(IOERR .GT. 0) GOTO 999 
OPEN (3 , FILE-LIST1 , IOSTAT=IOERR) 
IF(IOERR .GT. 0) GOTO 999 
OPEN(4 , FILE-LIST2 , IOSTAT=IOERR) 
IF(IOERR .GT. 0) GOTO 999 


C 

C STARTING POINT OF FM4 ANALYSIS LOOP FOR EACH DATA POINT 
C 


DO 111, J-l.K 

PRINT *, 'PROCESSING DATA POINT NO. : ' ,J 
READ(1 , 100) TIME(J) 

100 FORMAT(F10 . 2) 

READ(3 , 300) AMPLF , PHAS EF , ATIME 
300 FORMAT(A22,A22,A22) 

READ ( 4 , 400) FTIME , RTIME , NTIME , STIME 
400 FORMAT ( Al 8 , Al 8 , Al 8 , Al 8 ) 


C 

C OPEN LOOP DEPENDANT FILES 
C 


OPEN ( 5 , FI LE=- AMPLF , IOSTAT=IOERR) 

IF(IOERR .GT. 0) GOTO 999 

OPEN ( 6 , FILE-PHASEF , IOSTAT-IOERR) f 

IF(IOERR .GT. 0) GOTO 999 

OPEN (7 , FILE-ATIME , IOSTAT=IOERR) 

IF(IOERR .GT. 0) GOTO 999 

OPEN ( 8 , FILE- FTIME , IOSTAT-IOERR) 

IF(IOERR .GT. 0) GOTO 999 

OPEN(9 , FILE=RTIME , IOSTAT-IOERR) 

IF(IOERR .GT. 0) GOTO 999 

OPEN (10 , FILE-NTIME , IOSTAT-IOERR) 

IF( IOERR . GT . 0) GOTO 999 

OPEN( 11 , FILE-STIME , IOSTAT-IOERR) 


74 


APPB 


0111 


IF(I0ERR .GT. 0) GOTO 999 


0112 

c 



0113 

c 

READ IN PHASE AND FREQUENCY FILES 


0114 

c 



0115 


DO 222, 1-1,401 


0116 


READ (5 , 500) FREQ ( I ), FMAGRMS ( I ) 


0117 


500 F0RMAT(4X, F12 . 3 , 20X, F16 . 8) 


0118 


FMAG(I) - FMAGRMS ( I ) * (2. **.5) 


0119 


READ(6 , 600) PHAS(I) 


0120 


600 FORMAT (20X.F12. 3) 


0121 


222 CONTINUE 


0122 

c 



0123 

c 

PERFORM IFFT ON ORIGINAL SIGNAL 


0124 

c 



0125 


CALL IFOUR(FREQ , FMAG , PHAS ,TO,YO) 


0126 

c 



0127 

c 

PRINT TO FILE, 2 CYCLES OF ORIGINAL TIME 

SIGNAL, 

0128 

c 



0129 


NCYC2 - NCYC + NCYC 


0130 


DO 333, I-1.NCYC2 


0131 


WRITE(7 ,700) TO(I) ,YO(I) 


0132 


700 FORMAT (F16 . 8 , 4X, F16 .8) 


0133 


333 CONTINUE 


0134 

c 



0135 

c 

PERFORM FILTERING USING FILTERS A, B, C 


0136 

c 



0137 


IF(NN .LT. 1) GOTO 234 


0138 


CALL FILTER (FREQ , FMAG , FFMAG , AFILTS , AFILTE , BFILTS , BFILTE , 

0139 


*CFILTS , CFILTE) 


0140 

c 



0141 

c 

PERFORM IFFT ON FILTERED SIGNAL 


0142 

c 



0143 


CALL IFOUR(FREQ, FFMAG, PHAS, TO, FYO) 


0144 


234 CONTINUE 


0145 


IF(NN .LT. 1) THEN 


0146 


DO 345, 1-1,401 


0147 


FFMAG ( I ) - FMAG ( I ) 


0148 


345 CONTINUE 


0149 


DO 456, I-1.NCYC2 


0150 


FYO(I) - Y0(I) 


0151 


456 CONTINUE 


0152 


END IF 


0153 

c 



0154 

c 

PERFORM FILTERING TO FM4 SPECIFICATIONS 


0155 

c 



0156 


CALL FFM4(FREQ, FFMAG, PRIM, SPD.RMAG) 


0157 

c 



0158 

c 

PERFORM IFFT ON FM4 FILTERED SIGNAL 


0159 

c 

(USES ABC FILTERED SIGNAL AS BASIS) 


0160 

c 



0161 


CALL IFOUR(FREQ,RMAG ,PHAS .TO.RFYO) 


0162 

c 



0163 

c 

DETERMINE DIFFERENCE SIGNAL ( ORIGINAL - 

REGULAR COMPONENTS) 

0164 

c 



0165 


DO 444, 1=1, NCYC 
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0166 



DYO(I) = FYO(I) - RFYO(I) 

0167 


444 

CONTINUE 

0168 

c 



0169 

c 

DETERMINE SQUARED DIFFERENCE SIGNAL 

0170 

c 



0171 



DO 555, I-l.NCYC 

0172 



DSYO(I) - DYO(I)**2. 

0173 


555 

CONTINUE 

0174 

c 



0175 

c 

DETERMINE AVERAGE VALUE OF DIFFERENCE SIGNAL 

0176 

c 



0177 



DSUM - 0.0 

0178 



DO 666, I-l.NCYC 

0179 



DSUM - DSUM + DYO(I) 

0180 


666 

CONTINUE 

0181 



DBAR - DSUM/FLOAT (NCYC) 

0182 

c 



0183 

c 

DETERMINE STANDARD DEVIATION OF NET SIGNAL 

0184 

c 



0185 



CALL STDEV ( DYO , NCYC , DBAR , SDEV ) 

0186 

c 



0187 

c 

DETERMINE NORMALIZED KURTOSIS OF NET SIGNAL 

0188 

c 



0189 



CALL KURTOS (DYO, NCYC, DBAR, RKURT) 

0190 

c 



0191 

c 

WRITE RESULTS TO FILES 

0192 

c 



0193 



DO 777, I— 1 , NCYC 

0194 



WRITE(8 ,700) TO(I) ,FY0(I) 

0195 


777 

CONTINUE 

0196 

c 



0197 



DO 778, I— 1 , NCYC 

0198 



WRITE(9 , 700) TO(I) ,RFYO(I) 

0199 


778 

CONTINUE 

0200 

c 



0201 



DO 779, I-l.NCYC 

0202 



WRITE(10 , 700) TO ( I ) , DYO ( I ) 

0203 


779 

CONTINUE 

0204 

c 



0205 



DO 780, I-l.NCYC 

0206 



WRITE(11 , 700) T0(I) , DSYO(I) 

0207 


780 

CONTINUE 

0208 

c 



0209 



WRITE(2 , 200) TIME(J) , SDEV, RKURT 

0210 


200 

FORMAT ( FI 0 . 2 , F12 . 6 , F12 . 6) 

0211 

c 



0212 

c 

CLOSE LOOP DEPENDANT FILES 

0213 

G 



0214 



CLOSE ( 5 , IOSTAT-IOERR) 

0215 



IF( IOERR .GT. 0) GOTO 999 

0216 



CLOSE (6 , IOSTAT-IOERR) 

0217 



IF( IOERR .GT. 0) GOTO 999 

0218 



CLOSE ( 7 , IOSTAT=IOERR) 

0219 



IF( IOERR .GT. 0) GOTO 999 

0220 



CLOSE (8 , IOSTAT=IOERR) 
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0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 
0261 
0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 


IF(IOERR .GT. 0) GOTO 999 
CLOSE (9, IOSTAT-IOERR) 

IF(IOERR .GT. 0) GOTO 999 
CLOSE (10 , IOSTAT=IOERR) 

IF(IOERR .GT. 0) GOTO 999 
CL0SE(11 , IOSTAT-IOERR) 

IF(IOERR .GT. 0) GOTO 999 
C 

C CLOSE LOOP, MOVE TO NEXT DATA POINT 
C 

111 CONTINUE 
C 
C 

PRINT *, 'FM4 VALUES STORED IN FILE : \FM4F 
STOP 

999 PRINT *, 'ERROR NO. ' , IOERR 
END 
C 
C 

C END OF MAIN PROGRAM 

C 

C 

Q&*****ic***************icicjcici:****&&&&&&*&*ic*&*&**&*&******X*****'k'k*'k‘k*C 


C C 
C ROUTINE: I FOUR C 
C C 
C THIS SUBROUTINE PERFORMS AN INVERSE DISCRETE FOURIER C 
C TRANSFORM ON A DFT TRACE USING MAGNITUDE AND PHASE DATA. C 
C C 
C C 


c 

SUBROUTINE I FOUR (F , FMAG , PHAS , T , Y) 

C 

DIMENSION T(1024) ,Y(1024) 

DIMENSION F(401) ,DFTRE(401) ,DFTIM(401) ,FMAG(401) ,PHAS(401) 
DIMENSION FQ(20) 

COMPLEX X(1024) 

C 

C.. DEFINE DFT SETTINGS 
C 

NPTS - 1024 

NPOW = 10 

NFREQ - 401 

FSPAN - F(401) - F(l) 

PERIOD - FLOAT ( (NFREQ- 1) ) / FSPAN 
BW - F(2) - F(l) 

DELTA = PERIOD / FLOAT ( (NPTS-1) ) 

PRINT *, PERIOD 
PRINT *, BW 
C 

C INITIALIZE TIME, AND DFT ARRAYS 
C 

DO 100 1=1, NPTS 

T(I) = DELTA * FLOAT ( 1-1 ) 
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0276 



X(I) - ( 0..0. ) 

0277 

100 

CONTINUE 

0278 

c 



0279 

c 

CONVERTING PHASE AND MAGNITUDE TO REAL AND IMAGINARY FFT VALUES 

0280 

c 



0281 



PI - 3.1415927 

0282 



DO 300, I-l.NFREQ 

0283 



RAD - PHAS(I)*PI/180. 

0284 



DFTRE(I) - FMAG(I) * COS (RAD) 

0285 



DFTIM(I) - FMAG(I) * SIN(RAD) 

0286 


300 

CONTINUE 

0287 

c 



0288 

c 

..PERFORM INVERSE DFT ON REAL PART 

0289 

c 



0290 



I FLAG - 1 

0291 



DO 2200 I-l.NFREQ 

0292 



X(I) - CMPLX( DFTRE ( I ) , 0 . ) 

0293 

2200 

CONTINUE 

0294 



CALL FFT( I FLAG , NPOW , X ) 

0295 



DO 2300 I-l.NPTS 

0296 



Y(I) - REAL( X(I) ) 

0297 

2300 

CONTINUE 

0298 

C 



0299 

C 

. . PERFORM INVERSE DFT ON IMAGINARY PART 

0300 

C 



0301 



DO 2400 I-l.NFREQ 

0302 



X(I) - CMPLX( DFTIM(I) , 0 . ) 

0303 

2400 

CONTINUE 

0304 



DO 2500 I=NFREQ+1 , NPTS 

0305 



X(I) - ( 0..0. ) 

0306 

2500 

CONTINUE 

0307 



CALL FFT( I FLAG, NPOW, X ) 

0308 



DO 2600 I— 1 , NPTS 

0309 



Y(I) - Y(I) - AIMAG( X(I) ) 

0310 

2600 

CONTINUE 

0311 

c 



0312 



RETURN 

0313 



END 

0314 

c 



0315 

c 



0316 


0317 

c 


1 

0318 

c 


ROUTINE : FFT 1 

0319 

c 


i 

0320 

Q**********************************************' **********************{ 

0321 

c 



0322 



SUBROUTINE FFT( IFLAG.NPOW.X ) 

0323 



DIMENSION X(l) ,CS(2) ,MSK(13) 

0324 



COMPLEX X , CXC S , HOLD , XA 

0325 



EQUIVALENCE (CXCS,CS) 

0326 

c 



0327 



PI - 3.141 592 7 

0328 



NMAX - 2**NPOW 

0329 



ZZ - 2. * PI * FLOAT (I FLAG) / FLOAT (NMAX) 

0330 



MSK(l) - NMAX / 2 
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0331 

0332 

0333 

0334 

0335 

0336 

0337 

0338 

0339 

0340 

0341 

0342 

0343 

0344 

0345 

0346 

0347 

0348 

0349 

0350 

0351 

0352 

0353 

0354 

0355 

0356 

0357 

0358 

0359 

0360 

0361 

0362 

0363 

0364 

0365 

0366 

0367 

0368 

0369 

0370 

0371 

0372 

0373 

0374 

0375 

0376 

0377 

0378 

0379 

0380 

0381 

0382 

0383 

0384 

0385 


DO 15 1=2 , NPOW 

MSK(I) - MSK(I-l) / 2 
15 CONTINUE 

NN - NMAX 
MM - 2 
C 

C . . LOOP OVER NPOW LAYERS 
C 

DO 45 LAYER-1, NPOW 
NN - NN / 2 
NW = 0 

DO 40 1=1, MM, 2 
II - NN * I 
C 

C..CXCS - CEXP( 2*PI*NW*FLOAT(IFLAG) /FLOAT (NMAX) ) 

C 

W - FLOAT (NW) * ZZ 
CS(1) - COS( W ) 

CS (2) - SIN ( W ) 

C 

C . . COMPUTE ELEMENTS FOR BOTH HALVES OF EACH BLOCK 
C 

DO 20 J-l.NN 

II - II + 1 
IJ - II - NN 
XA = CXCS * X(1I) 

X(II) = X(IJ) - XA 
X(IJ) = X(IJ) + XA 
20 CONTINUE 

C 

C..BUMP UP SERIES BY 2, COMPUTE REVERSE ADDRESS 
C 

DO 25 LOC-2.NPOW 

LL = NW - MSK(LOC) 

IF( LL ) 30,35,25 

25 NW - LL 

30 NW - MSK(LOC) + NW 

GOTO 40 

35 NW - MSK(L0C+1) 

40 CONTINUE 

MM - MM * 2 
45 CONTINUE 

C 

C. .DO FINAL REARRANGEMENT, ALSO MULTIPLY BY DELTA. 

C 

NW - 0 

DO 80 1=1, NMAX 
NW1 - NW + 1 
HOLD - X(NW1) 

IF( NW1-I ) 60,55,50 

50 X(NW1) - X(I) 

55 X(I) = HOLD 

C 

C. .BUMP UP SERIES BY 1, AND COMPUTE REVERSE ADDRESS. 
C 
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0386 

60 

DO 65 LOC-l.NPOW 

0387 


LL = NW - MSK(LOC) 

0388 


IF( LL ) 70,75,65 

0389 

65 

NW - LL 

0390 

70 

NW = MSK(LOC) + NW 

0391 


GOTO 80 

0392 

75 

NW = MSK(LOC+l) 

0393 

80 

CONTINUE 

0394 


RETURN 

0395 


END 

0396 

C 


0397 

C 



0398 c*******^**********^************************^***^ r **'*'^ r, ® r ***************C 


0399 C C 

0400 C ROUTINE: FILTER C 

0401 C C 

0402 C THIS SUBROUTINE USES FILTERS A, B, AND C DEFINED IN THE MAIN C 

0403 C PROGRAM TO FILTER OUT ANY NON -GEAR -RELATED SIGNALS FROM THE C 

0404 C ANALYSIS. (BEARING NOISE, SLAVE GEAR NOISE, Etc.) C 

0405 C C 

0406 C C 


0407 

0408 C 

0409 SUBROUTINE FILTER (FREQ, FMAG , FFMAG , AFS , AFE.BFS , BFE , CFS , CFE) 

0410 C 

0411 DIMENSION FREQ(401) , FMAG(401) , FFMAG(401) 

0412 C 

0413 DO 100, 1=1,401 

0414 FFMAG ( I ) = FMAG(I) 


0415 

I F ( ( FREQ ( I ) .GE. 

AFS) 

.AND. 

(FREQ(I) 

.LE. AFE)) 

THEN 

0416 

0417 

0418 

FFMAG ( I ) - 0.0 
END IF 

IF( (FREQ(I) .GE. 

BFS) 

.AND, 

(FREQ(I) 

.LE. BFE)) 

THEN 

0419 

0420 

0421 

FFMAG ( I ) = 0.0 
END IF 

IF( (FREQ(I) .GE. 

CFS) 

.AND. 

(FREQ(I) 

.LE. CFE)) 

THEN 


0422 FFMAG ( I ) =0.0 

0423 END IF 

0424 100 CONTINUE 

0425 RETURN 

0426 END 

0427 C 

0428 C 

0429 c^^*^^*^****^^^*^*^^***^^*^^^^^*^^^^^^^***^**^^^^^*^ 


0430 C C 

0431 C ROUTINE: FFM4 C 

0432 C C 

0433 C THIS SUBROUTINE FILTERS OUT ALL AMPLITUDES EXCEPT: PRIMARY C 


0434 C MESH FREQUENCY AND FIRST HARMONIC AND THEIR FIRST CORRESPONDING C 

0435 C SIDEBANDS, AND PRIMARY SHAFT FREQ. (PER STEWART'S 1977 PAPER) C 

0436 C 

0437 C 

0438 C******************************************************************* 

0439 C 

0440 SUBROUTINE FFM4 (FREQ , FFMAG , PRIM , SPD ,RMAG) 
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0441 

c 



0442 



DIMENSION FREQ(401) ,FFMAG(401) ,RMAG(401) 

0443 

c 



0444 



HARM - PRIM * 2. 

0445 



LPRIM - PRIM - 1.5 * SPD 

0446 



HPRIM - PRIM + 1.5 * SPD 

0447 



LHARM - HARM - 1.5 * SPD 

0448 



HHARM - HARM + 1.5 * SPD 

0449 



LSHAFT = 0.5 * SPD 

0450 



HSHAFT - 1.5 * SPD 

0451 

c 



0452 



DO 100, 1=1,401 

0453 



RMAG(I) = 0.0 

0454 



IF( (FREQ(I) .GT. LPRIM) .AND. ( FREQ ( I ) .LT. HPRIM)) THEN 

0455 



RMAG(I) - FFMAG(I) 

0456 



END IF 

0457 



IF( (FREQ(I) .GT. LHARM) .AND. ( FREQ ( I ) .LT. HHARM)) THEN 

0458 



RMAG(I) = FFMAG(I) 

0459 



END IF 

0460 



IF( (FREQ(I) .GT. LSHAFT) .AND. (FREQ(I) .LT. HSHAFT)) THEN 

0461 



RMAG(I) = FFMAG(I) 

0462 



END IF 

0463 


100 

CONTINUE 

0464 



RETURN 

0465 



END 

0466 

c 



0467 

c 



0468 

C******************************************************************** c 

0469 

c 


c 

0470 

c 


ROUTINE : STDEV c 

0471 

c 


C 

0472 

c 


THIS SUBROUTINE COMPUTES THE STANDARD DEVIATION OF THE C 

0473 

c 


DIFFERENCE TIME SIGNAL (ORIGINAL -- REGULAR COMPONENTS). C 

0474 

c 


C 

0475 

c 


C 

0476 


0477 

c 



0478 



SUBROUTINE STDEV (DYO , N , DBAR , SD) 

0479 

c 



0480 



DIMENSION DYO(220) 

0481 



ASUM =0.0 

0482 



DO 100, 1=1, N 

0483 



ASUM = (DYO(I) - DBAR)**2. + ASUM 

0484 


100 

CONTINUE 

0485 



SD = (ASUM/FL0AT(N))**.5 

0486 



RETURN 

0487 



END 

0488 

c 



0489 

c 



0490 

£********************************************************************0 

0491 

c 


C 

0492 

c 


ROUTINE: KURTOS C 

0493 

c 


C 

0494 

c 


THIS SUBROUTINE COMPUTES THE NORMALIZED KURTOSIS OF THE C 

0495 

c 


DIFFERENCE TIME SIGNAL (ORIGINAL - REULAR COMPONENTS) C 
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0496 

c 





0497 

c 




c 

0498 


0499 

c 





0500 



SUBROUTINE KURTOS ( DYO , N , DBAR , RKT ) 


0501 

c 





0502 



DIMENSION DYO(220) 


0503 

c 





0504 



TSUM - 0.0 



0505 



DO 100, I— 1 , N 



0506 



TSUM = TSUM + 

(DYO(I) - DBAR)**4 . 


0507 


100 

CONTINUE 



0508 

c 





0509 



BSUM - 0.0 



0510 



DO 200, 1-1 ,N 



0511 



BSUM - BSUM + 

( DYO ( I ) - DBAR)**2 . 


0512 


200 

CONTINUE 



0513 

c 





0514 



RKT - ( FLOAT (N) * TSUM) /BSUM** 2 . 


0515 



RETURN 



0516 



END 
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0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 


C********************************************************************C 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

PROGRAM: HILB. FOR REVISION: A C 

REVISION DATE: 05-12-89 C 
PROGRAMMER: J.J. Zakrajsek DATE CREATED: 05-03-89 C 

THIS PROGRAM USES THE HILBERT TRANSFORM THEORY TO CONVERT 
THE TIME SIGNAL TO A COMPLEX FUNCTION TO DEMODULATE THE 
SIDEBANDS FROM THE CARRIER FREQUENCY. THIS METHOD IS BASED 
ON THE 1986 PAPER BY P. D. McFADDEN, AND THE 1988 PAPER BY 
D. A. WALLACE AND M. K. DARLOW. 

REV A: ADD VARIABLE RAMP FREQUENCY ANALYTIC PHASE VALUE 
CORRECTION FUNCTION. 




c 

DIMENSION TIME(50) ,AMAGN(220) ,APHAS(220) ,PPHAS(220) ,CPHAS(220) 
DIMENSION FREQ(401) ,FMAGRMS(401) ,FMAG(401) ,PHAS(401) ,FFMAG(401) 
DIMENSION T0(1024) ,YACT(1024) ,YHIL(1024) , POINT(60) 

CHARACTER*22 LIST1 , LIST3 , AMPLF, PHASEF 
CHARACTER*22 AAMPLF, APHASF 
COMPLEX X(1024) 

C 

C ENTER GENERAL PARAMETERS AND FILE NAMES 
C 

PRINT *, 'ENTER NUMBER OF DATA POINTS (INTEGER)' 

READ(*,*) K 

PRINT *, 'ENTER CARRIER FREQUENCY (PRIMARY OR FIRST HARMONIC)' 
READ ( * , * ) PRIM 

PRINT *, 'ENTER ROTATIONAL SPEED IN REV/SEC (Hz)' 

READ(*,*) SPD 

PRINT *, 'ENTER NUMBER OF TIME DATA POINTS PER CYCLE' 

PRINT *, '( INTEGER, 220 MAXIMUM )' 

READ (*,*) NCYC 
C 

C DEFINITION OF FREQUENCY BAND TO BE USED IN ANALYSIS 
C 

PRINT *, 'THIS ANALYSIS IS BASED ON A BAND OF FREQUENCY DATA' 
PRINT *, 'SURROUNDING THE MESH FREQUENCY CHOSEN ABOVE. THE ' 
PRINT *, 'NUMBER OF SIDEBANDS ON EACH SIDE OF THE FREQUENCY' 
PRINT *, 'OF INTEREST ARE CHOSEN BY THE USER. DATA OUTSIDE THIS' 
PRINT *, 'BAND, ALONG WITH THE FREQUENCY CHOSEN ABOVE, ARE' 

PRINT *, 'FILTERED OUT OF THE ANALYSIS.' 

PRINT *, ' 

C 

PRINT *, 'ENTER NUMBER OF SIDEBANDS TO BE PASSED BELOW' 

PRINT *, 'CHOSEN FREQUENCY ( REAL NUMBER XX.)' 

READ(* ,*) FILTS 

PRINT *, 'ENTER NUMBER OF SIDEBANDS TO BE PASSED ABOVE' 

PRINT *, 'CHOSEN FREQUENCY (REAL NUMBER XX.)' 

READ ( * , * ) FILTE 
C 

C INPUT PRIMARY FILE NAMES FILES 
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APPC 


0056 

0057 

0058 

0059 

0060 
0061 
0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 
0081 
0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 

0099 

0100 
0101 
0102 

0103 

0104 

0105 

0106 

0107 

0108 

0109 

0110 


/ 


c 

PRINT *, 'ENTER TXXLIST1.FIL INPUT FILES FILENAME' 

READ(* , ' (A22) ' ) LIST1 

PRINT *, 'ENTER TXXLIST3 . FIL OUTPUT FILES FILENAME' 

READ(* , ' (A22) ' ) LIST3 
C 

C OPEN PRIMARY FILES 
C 

OPEN (3 , FILE-LIST1 , IOSTAT-IOERR) 

IF(IOERR .GT. 0) GOTO 999 
OPEN (4 , FILE-LIST3 , IOSTAT=IOERR) 

IF(IOERR .GT. 0) GOTO 999 
C 

C STARTING POINT OF ANALYSIS LOOP FOR EACH DATA POINT 
C 

DO 111, J-l.K 

PRINT *, 'PROCESSING DATA POINT NO. : ' ,J 
READ (3, 300) AMPLF, PHASEF 
300 FORMAT(A22 , A22) 

READ (4, 400) AAMPLF, APHASF 
400 FORMAT (A2 2 , A22) 

C 

C OPEN LOOP DEPENDANT FILES 
C 

OPEN (5 , FILE=AMPLF , IOSTAT-IOERR) 

IF(IOERR .GT. 0) GOTO 999 
OPEN(6 , FILE-PHASEF, IOSTAT-IOERR) 

IF(IOERR .GT. 0) GOTO 999 

OPEN ( 8 , FILE-AAMPLF , IOSTAT=IOERR) 

IF(IOERR . GT. 0) GOTO 999 

OPEN ( 9 , FILE-APHASF , IOSTAT=IOERR) 

IF(IOERR .GT. 0) GOTO 999 
C 

C READ IN PHASE AND FREQUENCY FILES 
C 

DO 222, 1-1,401 

READ (5 , 500) FREQ ( I ) , FMAGRMS ( I ) 

500 FORMAT (4X, F12 . 3 , 20X, F16 .8) 

FMAG(I) - FMAGRMS ( I ) * (2. **.5) 

READ (6 , 600) PHAS(I) 

600 FORMAT (20X.F12. 3) 

222 CONTINUE 
C 

C PERFORM BANDPASS FILTERING 
C 

CALL FILTER( FREQ , FMAG , FFMAG , FILTS , FILTE , PRIM , SPD) 

C 

C PERFORM IFFT ON FILTERED FREQUENCY COMPONENTS OF ACTUAL SIGNAL 
C 

IFL - 1 

CALL I FOUR ( FREQ , FFMAG , PHAS , TO , YACT , I FL) 

C 

C OBTAIN HILBERT TRANSFORM OF ACTUAL SIGNAL 
C 

C PERFORM IFFT ON FREQUENCY COMPONENTS OF HILBERT TRANSFORM OF SIGNAL 
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0111 C 

0112 IFL - 2 

0113 CALL IF0UR(FREQ, FFMAG , PHAS ,T0,YHIL, IFL) 

0114 C 

0115 C DETERMINATION OF ANALYTIC FUNCTION MAGNITUDE AND PHASE 

0116 C 

0117 PI - 3.1415927 

0118 DO 444, I-l.NCYC 

0119 AMAGN(I) - (YACT(I)**2 . + YHIL(I)**2 .)**(. 5) 

0120 APHAS(I) - (ATAN2 (YHIL(I) , YACT(I) ) )*(180 . /PI) 

0121 444 CONTINUE 

0122 C 

0123 C CORRECT PHASE VALUES FOR NOMINAL GEAR MESH FREQUENCY 

0124 C 

0125 C ADD 180 DEGREES TO ALL PHASE VALUES 

0126 C 

0127 DO 210, I-l.NCYC 

0128 PPHAS(I) - APHAS(I) + 180. 

0129 210 CONTINUE 

0130 C 

0131 C DETERMINE FREQUENCY OF RAMP FUNCTION 

0132 C 

0133 NCT - 0 

0134 LIM - NCYC - 1 

0135 DO 211, I-l.LIM 

0136 N - I + 1 

0137 XCHK - PPHAS(I) - PPHAS(N) 

0138 IF (XCHK .GT. 200.) THEN 

0139 NCT - NCT + 1 

0140 POINT (NCT) = (TO(I) + T0(N))/2. 

0141 END IF 

0142 211 CONTINUE 

0143 SFREQ - FLOAT (NCT- 1)/(P0INT(NCT) - POINT(l)) 

0144 PRINT *, SFREQ 

0145 C 

0146 C DETERMINE INITIAL OFFSET OF RAMP FUNCTION 

0147 C 

0148 I - 0 

0149 212 CONTINUE 

0150 I - I + 1 

0151 N - I + 1 

0152 XCHK - P PHAS (I) - PPHAS(N) 

0153 IF (XCHK .GT. 200.) GOTO 213 

0154 GOTO 212 

0155 213 PERIOD - 1.0/SFREQ 

0156 OFFSET - PERIOD - (TO(I) + TO(N))/2.0 

0157 PRINT *, OFFSET 

0158 C 

0159 C CORRECT PHASE ANGLE VARIATION FOR NOMINAL TOOTH MESH FREQUENCY 

0160 C 

0161 DO 214, I-l.NCYC 

0162 I F ( I .EQ. 1) THEN 

0163 FC - 1.0 

0164 END IF 

0165 IF(I .EQ. 1) GOTO 379 
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0166 



N - I - 1 

0167 



XCHK - PPHAS(N) - PPHAS(I) 

0168 



IF(XCHK .GT. 200.) THEN 

0169 



M - I - 1 

0170 


377 

M - M + 1 

0171 



MN - M + 1 

0172 



CCHK - PPHAS(M) - PPHAS(MN) 

0173 



IF(CCHK .GT. 200.) GOTO 378 

0174 



GOTO 377 

0175 


378 

LL - I - 1 

0176 



ULIM - (TO(MN) + TO(M) )/2 . 0 

0177 



BLIM - (TO(I) + TO(LL) )/2 . 0 

0178 

c 


PRINT *, ULIM 

0179 

c 


PRINT *, BLIM 

0180 



CFREQ - 1.0/ (ULIM - BLIM) 

0181 



FC - CFREQ/SFREQ 

0182 



OFFSET - - 1 . 0*BLIM 

0183 



END IF 

0184 


379 

CONTINUE 

0185 



CPHAS(I) - PPHAS(I) - (360.)*SFREQ*FC*(TO(I) 

0186 


214 

CONTINUE 

0187 

c 



0188 

c 

WRITE RESULTS TO FILES 

0189 

c 



0190 



DO 777, I-l.NCYC 

0191 



WRITE(8 , 700) TO ( I ) , AMAGN ( I ) 

0192 


700 

FORMAT (F16 . 8 , 4X, F16 .8) 

0193 


111 

CONTINUE 

0194 

c 



0195 



DO 778, I-l.NCYC 

0196 



WRITE(9 , 700) TO(I) , CPHAS (I) 

0197 


778 

CONTINUE 

0198 

c 



0199 

c 

CLOSE LOOP DEPENDANT FILES 

0200 

c 



0201 



CLOSE (5 , IOSTAT-IOERR) 

0202 



IF(I0ERR .GT. 0) GOTO 999 

0203 



CLOSE (6 , IOSTAT-IOERR) 

0204 



IF(IOERR .GT. 0) GOTO 999 

0205 



CLOSE (8 , IOSTAT-IOERR) 

0206 



IF(IOERR .GT. 0) GOTO 999 

0207 



CLOSE(9 , IOSTAT-IOERR) 

0208 



IF(IOERR .GT. 0) GOTO 999 

0209 

c 



0210 

c 

CLOSE LOOP, HOVE TO NEXT DATA POINT 

0211 

c 



0212 


in 

CONTINUE 

0213 

c 



0214 

c 



0215 



PRINT *, 'VALUES STORED IN FILES ON DISK :A' 

0216 



STOP 

0217 


999 

PRINT *, 'ERROR NO. \IOERR 

0218 



END 

0219 

c 



0220 

c 




OFFSET) 
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0221 

0222 

0223 

0224 

0225 

0226 

0227 

0228 

0229 

0230 

0231 

0232 

0233 

0234 

0235 

0236 

0237 

0238 

0239 

0240 

0241 

0242 

0243 

0244 

0245 

0246 

0247 

0248 

0249 

0250 

0251 

0252 

0253 

0254 

0255 

0256 

0257 

0258 

0259 

0260 
0261 
0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 


C END OF MAIN PROGRAM 

C 

C 

c********************************************************************c 


c c 

C ROUTINE: I FOUR C 
C C 
C THIS SUBROUTINE PERFORMS AN INVERSE DISCRETE FOURIER C 
C TRANSFORM ON A DFT TRACE USING MAGNITUDE AND PHASE DATA. C 
C C 
C C 


C 

SUBROUTINE IFOUR(F , FMAG , PHAS , T , Y, IFL) 

C 

DIMENSION T(1024) ,Y(1024) 

DIMENSION F(401) ,DFTRE(401) ,DFTIM(401) ,FMAG(401) ,PHAS(401) 
DIMENSION FQ(20) 

COMPLEX X(1024) 

C 

C.. DEFINE DFT SETTINGS 
C 

NPTS - 1024 

NPOW - 10 

NFREQ - 401 

FSPAN - F(401) - F(l) 

PERIOD - FLOAT ( (NFREQ- 1) ) / FSPAN 
BW - F(2) - F(l) 

DELTA - PERIOD / FLOAT ( (NPTS-1) ) 

PRINT *, PERIOD 
PRINT *, BW 
C 

C INITIALIZE TIME, AND DFT ARRAYS 
C 

DO 100 I-l.NPTS 

T(I) - DELTA * FLOAT ( 1-1 ) 

X(I) - ( 0 . , 0 . ) 

100 CONTINUE 

C 

C CONVERTING PHASE AND MAGNITUDE TO REAL AND IMAGINARY FFT VALUES 
C 

PI - 3.1415927 
C 

C ACTUAL SIGNAL ALGORYTHM (IFL - 1) 

C 

IF(IFL .EQ. 1) THEN 
DO 300, 1-1, NFREQ 
RAD - PHAS (I)*PI/180 . 

DFTRE(I) - FMAG(I) * COS(RAD) 

DFTIM(I) - FMAG(I) * SIN(RAD) 

300 CONTINUE 
END IF 
C 

C HILBERT TRANSFORM ALGORYTHM (IFL - 2) 

C 
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0276 

0277 

0278 

0279 

0280 
0281 
0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 

0301 

0302 

0303 

0304 

0305 

0306 

0307 

0308 

0309 

0310 

0311 

0312 

0313 

0314 

0315 

0316 

0317 

0318 

0319 

0320 

0321 

0322 

0323 

0324 

0325 

0326 

0327 

0328 

0329 

0330 


IF(IFL .EQ. 2) THEN 
DO 303, I-l.NFREQ 
RAD - PHAS (I)*PI/180 . 

DFTRE(I) - FMAG(I)*SIN(RAD) 

DFTIM(I) - ( -1 . 0)*FMAG(I)*COS(RAD) 

303 CONTINUE 
END IF 
C 

C. .PERFORM INVERSE DFT ON REAL PART 
C 

I FLAG - 1 
DO 2200 I-l.NFREQ 

X(I) - CMPLX( DFTRE ( I ) , 0 . ) 

2200 CONTINUE 

CALL FFT( I FLAG ,NPOW,X ) 

DO 2300 I-l.NPTS 

Y(I) - REAL( X(I) ) 

2300 CONTINUE 

C 

C.. PERFORM INVERSE DFT ON IMAGINARY PART 
C 

DO 2400 I-l.NFREQ 

X(I) - CMPLX( DFTIM(I) , 0 . ) 

2400 CONTINUE 

DO 2500 I-NFREQ+1 , NPTS 
X(I) - ( 0..0. ) 

2500 CONTINUE 

CALL FFT( IFLAG,NPOW,X ) 

DO 2600 I-l.NPTS 

Y(I) - Y(I) - AIMAG( X(I) ) 

2600 CONTINUE 

C 

RETURN 

END 

C 

C 

C***************************************************************** , ' r '**C 

c C 

C ROUTINE: FFT 

C 

c 

SUBROUTINE FFT( IFLAG,NPOW,X ) 

DIMENSION X(l) ,CS(2) ,MSK(13) 

COMPLEX X, CXCS, HOLD ,XA 

EQUIVALENCE (CXCS.CS) 

C 

PI - 3.141 592 7 
NMAX - 2**NPOW 

ZZ - 2. * PI * FLOAT (I FLAG) / FLOAT (NMAX) 

MSK(l) - NMAX / 2 
DO 15 I-2.NPOW 

MSK(I) - MSK(I-l) / 2 
15 CONTINUE 

NN - NMAX 
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0331 

MM 

- 2 


0332 

C 



0333 

C . . LOOP OVER NPOW LAYERS 


0334 

c 



0335 

DO 

45 LAYER-1, NPOW 


0336 


NN - NN / 2 


0337 


NW - 0 


0338 


DO 40 1-1 ,MM, 2 


0339 


II - NN * I 


0340 

C 



0341 

C. .CXCS - CEXP( 2*PI*NW*FLOAT ( I FLAG) /FLOAT (NMAX) ) 

0342 

c 



0343 


W - FLOAT (NW) * ZZ 


0344 


CS(1) - C0S( W ) 


0345 


CS (2) - SIN( W ) 


0346 

c 



0347 

C . , COMPUTE ELEMENTS FOR BOTH HALVES OF EACH 

BLOCK 

0348 

C 



0349 


DO 20 J-l.NN 


0350 


II - II + 1 


0351 


IJ - II - NN 


0352 


XA - CXCS * X(II) 


0353 


X(II) - X(IJ) - XA 


0354 


X(IJ) - X(IJ) + XA 


0355 

20 

CONTINUE 


0356 

C 



0357 

C. .BUMP UP SERIES BY 2, COMPUTE REVERSE ADDRESS 

0358 

C 



0359 


DO 25 LOC— 2 , NPOW 


0360 


LL - NW - MSK(LOC) 


0361 


IF( LL ) 30,35,25 


0362 

25 

NW - LL 


0363 

30 

NW - MSK(LOC) + NW 


0364 


GOTO 40 


0365 

35 

NW - MSK(LOC+l) 


0366 

40 

CONTINUE 


0367 


MM - MM * 2 


0368 

45 

CONTINUE 


0369 

C 



0370 

C..DO FINAL REARRANGEMENT, ALSO MULTIPLY BY 

DELTA . 

0371 

C 



0372 

NW 

- 0 


0373 

DO 

80 I-l.NMAX 


0374 


NW1 - NW + 1 


0375 


HOLD - X(NW1) 


0376 


IF( NW1-I ) 60,55,50 


0377 

50 

X(NW1) - X(I) 


0378 

55 

X(I) - HOLD 


0379 

C 



0380 

C . . BUMP UP SERIES BY 1 , AND COMPUTE REVERSE 

ADDRESS 

0381 

C 



0382 

60 

DO 65 L0C-1.NP0W 


0383 


LL - NW - MSK(LOC) 


0384 


IF( LL ) 70,75,65 


0385 

65 

NW - LL 
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0386 

70 

NW - MSK(LOC) + NW 

0387 


GOTO 80 

0388 

75 

NW - MSK(LOC+l) 

0389 

80 

CONTINUE 

0390 


RETURN 

0391 


END 

0392 

C 


0393 

C 



0394 C ******************************************************************** C 

0395 C 

0396 C ROUTINE: FILTER 

0397 C 

0398 C THIS SUBROUTINE FILTERS OUT ALL FREQUENCY COMPONENTS 

0399 C OUTSIDE THE FREQUENCY BAND CHOSEN IN THE MAIN PROGRAM. 

0400 C THE CARRIER FREQUENCY IS ALSO FILTERED OUT. 

0401 C 

0402 C 

0403 c******************************************************************* 1 ^ 

0404 C 

0405 SUBROUTINE FILTER ( FREQ , FMAG , FFMAG , FS , FE, PRIM, SPD) 

0406 C 

0407 DIMENSION FREQ(401) , FMAG(401) , FFMAG(401) 

0408 C 

0409 SBL - PRIM - (FS + .5)*SPD 

0410 PL - PRIM - ( . 5*SPD) 

0411 PU - PRIM + ( . 5*SPD) 

0412 SBU - PRIM + (FE + .5)*SPD 

0413 C 

0414 DO 100, 1-1,401 

0415 FFMAG ( I ) - 0.0 

0416 IF( (FREQ(I) .GE. SBL) .AND. (FREQ(I) .LE. PL)) THEN 

0417 FFMAG ( I ) - FMAG(I) 

0418 END IF 

0419 • IF( (FREQ(I) .GE. PU) .AND. ( FREQ ( I ) .LE. SBU)) THEN 

0420 FFMAG ( I ) - FMAG (I) 

0421 END IF 

0422 100 CONTINUE 

0423 RETURN 

0424 END 
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0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 


APPENDIX D - SOURCE CODE FOR PARAM. FOR 


c c 
c c 

C PROGRAM: PARAM. FOR CREATED: 05-05-89 REVISED: XX-XX-XXXX C 
C C 
C PROGRAMMER: J.J. ZAKRAJSEK C 
C C 
C THIS PROGRAM COMPUTES SEVERAL FAULT DETECTION PARAMETERS C 
C BASED ON THE A-B-C FILTERED SIGNAL CALCULATED IN THE PROGRAM C 
C TFM4 . FOR . THE BASIS OF THE PARAMETERS ARE SUMMARIZED BELOW. C 
C C 
C CREST FACTOR (CF) - PEAK LEVEL / TOTAL RMS LEVEL C 
C C 
C SIDEBAND LEVEL FACTOR (SBLF) - FIRST ORDER SIDEBANDS (PRIMARY) C 
C TOTAL RMS LEVEL C 
C C 
C WEAR ENERGY RATIO - RMS LEVEL OF NON-PERIODIC PART OF SIGNAL C 
C (WER) / RMS LEVEL OF PERIODIC PART OF SIGNAL C 
C C 
C C 


C 

DIMENSION TIME(50) ,FSIG(220) ,RSIG(220) ,DSIG(220) 

DIMENSION FREQ (500) , AMPL(500) 

CHARACTERS 2 FMOF, PARAMF , LI ST2 , LIST1 , AMPLF 
CHARACTER*18 FTIMEF , RTIMEF , DTIMEF 
C 

C ENTER INITIAL DATA AND FILE NAMES 
C 

PRINT *, 'ENTER NUMBER OF DATA POINTS' 

READ(*,*) K 

PRINT *, 'ENTER MESH FREQ (XXXX.)' 

READ(* , *) PRIM 

PRINT *, 'ENTER ROTATIONAL SPEED (REV/SEC)' 

READ(*,*) SPD 
C 

PRINT *, 'ENTER OUTPUT DATA FILE NAME (\TXXPARAM.DAT) ' 

READ (*, ' (A22) ' ) PARAMF 

OPEN(7 , FILE=PARAMF , IOSTAT=IOERR) 

IF ( IOERR .GT. 0) GOTO 999 

PRINT *, 'ENTER FILE NAMES ARRAY FILE (\TXXLIST1 . FIL) ' 

READ(* , ' (A22) ' ) LIST1 

OPEN (6 , FILE=LIST1 , IOSTAT^IOERR) 

IF(I0ERR .GT. 0) GOTO 999 

PRINT *, 'ENTER FILE NAMES ARRAY FILE (\TXXLIST2 . FIL) ' 

READ(*, ' (A22) ') LIST2 

OPEN (8 , FILE=LIST2 , IOSTAT=IOERR) 

IF(IOERR .GT. 0) GOTO 999 

PRINT *, 'ENTER PREVIOUS TFMO FILE (\TXXFMO.DAT)' 

READ(* , ’ (A22) ' ) FMOF 

OPEN(9 , FILE=FM0F, IOSTAT=IOERR) 

IF( IOERR .GT. 0) GOTO 999 

PRINT *, 'ENTER NUMBER OF TIME DATA POINTS PER CYCLE (INTEGER)' 
READ(*,*) NCYC 
C 
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0056 C START TIME LOOPING, PROCESSING DATA AT EACH TIME INTERVAL 

0057 C 

0058 DO 300, J-l.K 

0059 PRINT *, 'ACCESSING AND PROCESSING DATA FILE NO: ' ,J 

0060 READ(6 , 600) AMPLF 

0061 600 FORMAT (A20) 

0062 READ(8 , 800) FTIMEF , RTIMEF , DTIMEF 

0063 800 FORMAT(Al8 , Al8 , A18) 

0064 READ(9 , 900) TIME(J) 

0065 900 FORMAT (F8. 2) 

0066 OPEN (4 , FILE-FTIMEF , IOSTAT-IOERR) 

0067 IF (IOERR .GT. 0) GOTO 999 

0068 OPEN ( 3 , FILE-RTIMEF , IOSTAT-IOERR) 

0069 IF (IOERR .GT. 0) GOTO 999 

0070 OPEN (2 , FILE-DTIMEF, IOSTAT-IOERR) 

0071 IF (IOERR .GT. 0) GOTO 999 

0072 OPEN ( 5 , FILE-AMPLF, IOSTAT-IOERR) 

0073 IF (IOERR .GT. 0) GOTO 9?9 . 

0074 C 

0075 C READ IN AND FIND THE MAX VALUE OF THE FILTERED TIME RECORD 

0076 C 

0077 PKL - 0.0 

0078 DO 100, I-l.NCYC 

0079 READ(4 , 400) FSIG(I) 

0080 400 FORMAT (20X, FI 6. 8) 

0081 AVALU - ABS(FSIG(I) ) 

0082 IF( AVALU .GT. PKL ) THEN 

0083 PKL - AVALU 

0084 END IF 

0085 100 CONTINUE 

0086 C 

0087 C FIND FIRST ORDER SIDEBANDS LEVEL 

0088 C 

0089 DO 301, 1-1,401 

0090 READ(5 , 500) FREQ(I) ,AMPL(I) 

0091 500 FORMAT (4X, F12 . 3 , 20X, F16 . 8) 

0092 301 CONTINUE 

0093 C 

0094 CALL FOSB ( FREQ, AMPL, PRIM, S PD, SBL) 

0095 C 

0096 C FIND TOTAL RMS LEVEL OF SIGNAL 

0097 C 

0098 CALL STDEV(FSIG ,NCYC , FSD) 

0099 C 

0100 C FIND RMS LEVEL OF PERIODIC COMPONENTS 

0101 c 

0102 DO 120, I-l.NCYC 

0103 READ (3 , 400) RSIG(I) 

0104 120 CONTINUE 

0105 C 

0106 CALL STDEV (RSIG , NCYC , RSD) 

0107 C 

0108 C FIND RMS LEVEL OF DIFFERENCE SIGNAL (SIGNAL - PERIODIC COMPONENTS) 

0109 C 

0110 DO 140, I-l.NCYC 
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0111 READ(2 , 400) DS1G(I) 

0112 140 CONTINUE 

0113 C 

0114 CALL STDEV(DSIG ,NCYC ,DSD) 

0115 C 

0116 C COMPUTE CREST FACTOR 

0117 C 

0118 CF - PKL/FSD 

0119 C 

0120 C FIND SIDEBAND LEVEL FACTOR 

0121 C 

0122 SBLF - SBL/FSD 

0123 C 

0124 C FIND WEAR ENERGY RATIO 

0125 C 

0126 WER - DSD/RSD 

0127 C 

0128 C CLOSE LOOP DEPENDANT FILE, AND WRITE CALCULTED DATA TO FILE 

0129 C 

0130 CLOSE (4 , I0STAT=I0ERR) 

0131 IF (IOERR .GT. 0) GOTO 999 

0132 CLOSE (2 , IOSTAT—IOERR) 

0133 IF (IOERR .GT. 0) GOTO 999 

0134 CLOSE (3 , IOSTAT-IOERR) 

0135 IF (IOERR .GT. 0) GOTO 999 

0136 CLOSE(5 , IOSTAT-IOERR) 

0137 IF (IOERR .GT. 0) GOTO 999 

0138 WRITE (7,700) TIME(J) ,CF, SBLF, WER 

0139 700 FORMAT (F8 . 2 , F10 .4 , F10 . 4 , F10 .4) 

0140 300 CONTINUE 

0141 C 

0142 PRINT *, 'PARAMETER VALUES STORED IN FILE ; ' , PARAMF 

0143 STOP 

0144 999 PRINT*, 'ERROR NO. ', IOERR 

0145 END 

0146 C 

0147 C 

0148 C ****** END OF MAIN PROGRAM ******* 

0149 C 

0150 C 

0151 C********************************************************************C 

0152 C 

0153 C ROUTINE: FOSB (FIRST ORDER SIDEBANDS) 

0154 C 

0155 C********************************************************************C 

0156 C 

0157 SUBROUTINE FOSB ( FREQ, AMPL, PRIM, SPD.SBL) 

0158 C 

0159 DIMENSION FREQ (500) ,AMPL(500) 

0160 BSBS - PRIM - 1.5*SPD 

0161 BSBE - PRIM - . 5*SPD 

0162 USBS - PRIM + . 5*SPD 

0163 USBE - PRIM + 1 . 5*SPD 

0164 BSB - 0.0 

0165 USB - 0.0 


93 


non 



APPD 


0166 

0167 

0168 

0169 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0178 

0179 

0180 
0181 
0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 

0197 

0198 

0199 

0200 


DO 200, 1-1,401 

IF((FREQ(I) .GT.BSBS) .AND. (FREQ(I) . LT . BSBE) ) THEN 
IF(AMPL(I) .GT. BSB) BSB - AMPL(I) 

ENDIF 

IF( (FREQ(I) .GT.USBS) .AND. (FREQ(I) .LT.USBE) ) THEN 
IF(AMPL(I) .GT. USB) USB - AMPL(I) 

ENDIF 

200 CONTINUE 

SBL - (BSB + USB)/2 . 0 
RETURN 
END 
C 
C 

c ********************************************************************c 
c c 

C ROUTINE: STDEV (STANDARD DEVIATION OF TIME SIGNAL) C 

C C 

C********************************************************************C 

c 

SUBROUTINE STDEV (DYO,N, SD) 

C 

DIMENSION DYO(220) 

DSUM - 0.0 
DO 100, I-l.N 
DSUM - DSUM + DYO(I) 

100 CONTINUE 

DBAR - DSUM/FLOAT (N) 

C 

ASUM - 0.0 
DO 200, I-l.N 

ASUM - (DYO(I) - DBAR)**2 . + ASUM 
200 CONTINUE 

SD - ( ASUM/FLOAT (N) ) ** . 5 

RETURN 

END 
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A amplitude at mesh frequency and harmonics 

A(t) original signal after time synchronous averaging 

AN(f ) analytic signal 

D(t) difference signal 

d discrete value In the difference signal, D(t> 

mean value of the difference signal, D(t) 

F(m> discrete value In a frequency array 

f(n) discrete value In a time record 

F(m) discrete frequency array (discrete fourler transform of a time 
record) 

f(n) discrete time record 

F(w) frequency domain function (fourler transform of a time domain func- 
tion) 

f frequency (cycles/second) 

f(mln) minimum resolution frequency 

FOSL first order sidebands level 

Gxy(f) cross power spectrum 

Gxy*(f) complex conjugate of Gxy(f) 

Gyx(f) cross power spectrum 

Gxx(f) auto power spectrum of Input 

Gyy(f) auto power spectrum of output 

H(f) frequency response function 

Hyx(f) estimate of frequency response function 

h(t) unit Impulse response function 

K kurtosls, fourth moment of an array of values 

N number of data points In a discrete time record 

NK normalized kurtosls 
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NS number of samples 

Nt number of time averages 

PL peak level of signal 

PP peak to peak level of signal average 

R(t) regular meshing components of signal average 

RMS standard deviation of signal average 

RMSDS standard deviation of difference signal, D(t) 

RMSRC standard deviation of regular meshing components of signal, R(t) 

S/N signal to noise ratio 

T length of time record (seconds) 

t time (seconds) 

w frequency (radians/second) 

X(f) fourler transform of Input signal 

X*(f ) complex conjugate of X(f) 

Y(f) fourler transform of output signal 
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